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Abstract 


Geysering is a type of unstable and periodic boiling occurring during start-up. It causes flow 
oscillations which change the void fractions and reactivity. This makes the nuclear reactor 
difficult to control. Hence bench-marking of geysering is required. In the present code, a four 
equation Drift-Flux model has been used for better numerical stability. Each component has 
been represented one dimensionally with variable cross-sectional area. The equations are 
solved by a partially implicit finite difference technique that can use different time steps in 
different components. The components are discretised using staggered mesh arrangements. In 
the finite difference approximation explicit updating of velocity has been done while pressure 
calculation has been done implicitly. 

At onset of nucleate boiling (ONE), the bubble formation rate has been calculated from 
the evaporation heat consumed during growth and the latent heat of vaporization. The point at 
which the bubbles are first detached from the wall (this is OSV: onset of significant void) is 
determined by the sub cooling of the liquid which is dependent upon the thermal and 
hydrodynamic condition in the channel. The departure radius of the bubble at this point has 
been calculated using Unal’s semi-empirical model. The growth and condensation behavior of 
small bubbles are inertia controlled and that of large bubbles are governed by heat transfer at 
the phase interface. Jacob number has been used to decide whether the flow is inertia dominant 
or controlled by heat transfer. 

Inertia controlled condensation has been estimated by Hamit’s correlation. Heat 
transfer controlled condensation has been calculated using Hainoun et al’s model. The model 
incorporated in this work determines the void fraction in the sub cooled boiling regime. The 
simulation has been verified by comparing with the McMaster experiments on axial void 
distribution. 

Finally the code has been run and validated for the experimental Apsara Loop at 
Bhabha Atomic Research Centre. The pressure drop in the single phase leg has been 
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satisfactorily matched with the experimental data. The simulated void fraction distribution in 
the two phase leg has been attempted to match with experimental results. But there are still 
some anomalies while choosing the thermal constraint for the energy equation. The 
assumptions like thermal equilibrium of the two phases or saturated condition of one of the 
phases, while considering mixture energy equation, are very crucial to the final solution. 
Above all a proper boundary condition, proper convergence criteria in the pressure iteration 
loop and a judicial setting of datum value of void fraction (a^)is essential in getting a correct 

solution. The following important observations have been drawn from the present work: 

1) The code can simulate the small break LOCA satisfactorily. 

2) Inclusion of vapor generation and bubble condensation model was helpful in predicting the 
geysering phenomenon. The model has been validated with McMaster experimental test 
results on axial void distribution. 

3) The usefulness of the code has been studied for validating the Apsara Lx)op experimental 
test data for void fraction in the two phase leg at two locations, pressure drop in the single 
phase leg and the temperature in the loop at different locations. 
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Introduction 


Chapter 1 

Introduction 

1.1 Natural Circulation Loop 

In a natural circulation loop the circulating fluid removes the heat from a heat source and 
transports it to heat sink, the fluid circulation being established due to the buoyancy force. The 
buoyancy force is the driving force. Buoyancy force is result of thermally induced density 
differences in a gravitational field. The heat sink is located at a higher elevation than the heat 
source. Such loops find application in many engineering fields such as, gas turbine blade 
cooling, solar water heater, geothermal systems, electrical machine rotor cooling, transformer 
cooling, nuclear reactor core cooling etc. 

Present generation reactors like the Light Water Reactors (LWRs) and PHWR 
(Pressurized Heavy Water Reactors) utilize forced circulation loops for heat transfer from core 
to the steam generator. But this method has the following disadvantages: 

1 ) Usage of pump is costly 

2) If pump break down occurs then dissipation of fission heat is hampered. This results 
in tremendous accumulation of heat and may cause core melt down in the reactor. 
There are many examples of nuclear accidents due to core melt down like the Three 
Mile Island Unit-2 (1979) , Chernobyl (1986) and some less significant melting 
events like EBR-1 (National Reactor Testing Station, Idaho, USA, 1951), Windscale 
(now Sellafield, England, 1957), Lucens (Switzerland, 1969) and St. Laurent (France, 
1969). 

3) Back up safety measures have to be used. 
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Hence next generation of boiling water reactors should have natural circulation of 
primary coolant under normal operational situations. Figure 1.1 shows the schematic diagram 
of a two-phase natural circulation loop. The loop essentially consists of a tubular heated 
section, a riser, a steam drum or a separator, a condenser, a down comer, and upper- lower 
horizontal sections. The steam collected in separator or steam drum goes to the condenser, 
where it is condensed and then sent back to the separator. This loop differs from forced 
circulation loop only by the absence of pump. 

In a two-phase natural circulation loop since heat is removed by natural circulation, it 
is inherently safe. This method of cooling enhances passive safety. 

Many proposed advanced designs of nuclear power reactors use two-phase natural 
circulation of primary coolants. However a major problem with these loops is the occurrence 
of flow instabilities. How instabilities are undesirable in boiling, condensing and other two- 
phase flow systems for several reasons. 


1.2 Boussinesq Approximation 


In natural convection flows, the basic driving force arises from the temperature or 
concentration field. The temperature variation results in a difference in density. This in turn 
results in a buoyancy force due to the presence of the body force field. For a gravitational 
field, the body force F = pg, where, g is the force per unit mass of the fluid. It is the variation 
of p that gives rise to the flow and if this variation were to be neglected, no flow would 
result. The temperature field (energy equation) is linked with the flow and all of the equations, 
mass balance, momentum balance and energy balance, are coupled through the variation of the 
density p. Therefore, these equations have to be solved simultaneously to give the 
distributions, in space and time, of the velocity, pressure and temperature fields. Due to this 
added complexity in the analysis of flow, several simplifying assumptions and approximations 
are made in natural convection, to allow a more convenient procedure for obtaining a solution 
to the Navier-Stokes equation. 
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In the momentum equation, the local static pressure p may be broken down into two 
terms, one due to hydrostatic pressure, and the other due to the motion of the fluid, p^ . The 

former coupled with the body force acting on the fluid constitutes the driving mechanism for 
the flow. 

P = P„ + Pd (1-1) 


The pressure p^ is the hydrostatic pressure in the ambient medium. Therefore, for a 
gravitational field, it is given as, = P„8 ■> where g is the force due to gravity per unit 
mass of the fluid and is the density of the ambient fluid. Considering, for simplicity, an 
isothermal flat vertical surface at temperature Tj and immersed in an isothermal ambient 
medium r„. The driving buoyancy force for the vertical flow arises due to the difference 
between the body force and the force due to the hydrostatic pressure gradient in the ambient 
medium. This implies that the driving force per unit volume of the fluid isg(Pg - p). If this 
force acts over a vertical distance y and if the energy thus added is equated to the kinetic 

energy per unit volume, we have, ^ ^ = gy{p„ - p) 


Since the viscous forces have been neglected here and only the driving force has been 
considered, the above estimate of is the maximum value that may be expected in the flow. 

Considering an inviscid flow over the vertical surface, Bernoulli’s equation may be applied to 
estimate the order of magnitude of the pressure difference in the flow. If p^ is the pressure in 
the ambient medium and p is the pressure at a location in the boundary layer, the value of 
(Pa ~ p) be determined, knowing that the velocity in the ambient medium is zero. 

Therefore, 



and p„ - p « o\gy{p„ - p)] 


( 1 . 2 ) 


3 



Introduction 


y 



(1.3) 


Hence it is now necessary to determine the density difference in terms of the temperature and 
pressure in the flow. 

If the density of the fluid p is a function of the pressure and temperature, the density 
at a given point in the flow, p{p,T), may be written in terms of the density p„ in the 

ambient medium, as a double Taylor series expansion about the ambient conditions. 

This may be given as follows: 


P„ =P + 




(T.-T)+- 


l( d^p 




{T,-Tf+. 


+ 


\^Pj 
2 


(Pa-p)+:^hrT (Pa-pf + 


2 ! 


dp 




Now the coefficient of thermal volumetric expansion )3 is given by 

P[»T), 


(1.4) 


and. Pa- p~ o\gy{pa - p)] is the estimate of the maximum pressure difference. Therefore, 

Pi8^ 


P, - P = Pi3(r -r„ )+i^(r -r, f + . 


.^p]t 

2 


[?}’(p«-p)]+^f|-yl [gy(p.-p)f + 


2 ! 


dp 


JT 


+^[gy(p.-pXT,-T)h. 


(1.5) 


In the first series, the leading term is p/3 (T - ) 

The second term is ^-— (T-T ) 
2! 
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Hence the n'*’ term will be — - (7 - )" 

«! 

The series will be convergent when [j8(7 )]< 1 and if ^3(7 -7^ )]« 1, then only the first 
term of the series may be retained. 


If 





T 


« 1 , then the second series may be neglected. Similarly, 


dp 


8y 


L\ )T J 


and 


[/3(7 -7„ )] must both be much less than unity to neglect the third series. 

'(dp'^ 


Therefore, if |j8(7-7„)]«l and 
estimated as: p^- p = p^{r -7^ ) 


dp 


L'. JT 


gy 


«1 then the density difference may be 


General momentum equation for natural convection: 

DV 


=F-Vp+^vV+^v(v-\/) 


( 1 . 6 ) 


D 


where the operator — represents the substantial or particle derivative and given by 

D a ,, a, , a , a 
— = — + v — + v — + v — 

Dt dt ^dx "ay ^ dz 

V = Velocity vector = iV^ + JV^ + kV, 

F = Body force per unit volume = iF^ + jF^ + kF, 

d d d 

V = Vector Operator Del = i— + j — + ^ 

a^: dy dz 

pL = Coefficient of viscosity (dynamic viscosity) of the fluid 
The shear stresses and normal stresses are related as follows; 




av av 


dy dx 


= '^.r 


(1.7a) 
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(1.7b) 


f 9K "l 

1 dz ox 


(1.7c) 


Normal stress, = -p + 2{i — - + AV • V 

3x 


For an incompressible fluid, with constant density, V • V = 0 and therefore, 

o-^=-p + 2;t— ^ 


(1.9a) 


<7„ =-p + 2/r- 


(1.9b) 


a^=-p + 2p- 


(1.9c) 


A is known as the second coefficient of viscosity and for variable density and monatomic gas 
it is given by, A = -.j/i .Here , ^V^V + -^v(v.v) represents the viscous dissipation term 
culminating from shear stresses and normal stresses. 

Neglecting the viscous dissipation term, for a gravitational field, the body force F = pg 


Where g is the force per unit mass of the fluid (acceleration due to gravity). 

Local static pressure p = p„ + Pa 

p^ = Hydrostatic pressure in the ambient medium 
p^ = Dynamic component of pressure due to the flow of fluid 
For a gravitational field, = Pa8 


( 1 . 10 ) 


From the momentum equation: 


Or, p„ = p„gy 


F-Vp = pg-V(p„ +pj 


( 1 . 11 ) 
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= (P5-VpJ-Vp^ 

= (P5-P„^)-Vprf =(p-pJg-Vp^ (1.12) 

If g is downwards and y-direction upwards, as is the generally case for vertical buoyant 
flows, g = -jg , and hence 


F-Vp = (p„-p)s/-Vp^ (1.13) 

where j is the unit vector in y-direction and “g” is the magnitude of the gravitational force, 
per unit mass of the fluid (acceleration due to gravity). For g at an angle 6 with the y- 
direction, the term, therefore, becomes (p„ - p)g cos 6 in the y-direction and (p^ -p)gsin0 
in the x, or horizontal direction. 


pi ^ 

In Cartesian co-ordinates it becomes: - pg — — = {p„-p)g- 

dy dy 


= (p«-pk' 


d{p-Pn) 

dy 


Butp„-p = pP{T-Tj 

Hence, -pg -|P = pg/3(r-rj-fc^ 
ay ay 


(1.14) 

(1.15) 

(1.16) 


1.3 Classification of Flow Instabilities 


Flow instability phenomena can be classified on the basis of following definitions: 

A steady flow is one in which the system parameters are functions of the space variables only. 
Practically, however, they undergo small fluctuations due to turbulence, nucleation, or slug 
flow. These fluctuations play a role in triggering several instability phenomena. 

A flow is stable if, when momentarily disturbed, its new operating conditions tend 
asymptotically towards the initial ones. The various parameters that may affect the stability 
and can lead to instabilities are: 

Geometry-which includes channel length, channel diameter, inlet and outlet restrictions, single 
or multiple channels, 
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Operating conditions- such as pressure, inlet sub cooling, mass velocity, power input, forced or 
natural circulation. 

Boundary conditions- like axial heat flux distribution, pressure drop across the channels. 

Boure et al [1973] divided these instabilities into two classes: static instability and 
dynamic instability. 

Static instability is one where a small change in one of the independent variables leads to a 
large change in one of the dependent variables. A static instability can lead either to a different 
steady state condition or to a periodic behavior. Static instabilities can be described by steady 
state equilibrium also. Examples of static instability are listed below: 

Flow regime transition: Cyclic flow pattern transitions and flow pattern variation occurs. 
Boiling crisis: Bubbly flow has less void but higher pressure drop than that of annular flow 
characterized by the ineffective removal of heat from the heated surface, leads to wall 
temperature excursion and flow oscillation. 

The Ledinegg instability or flow excursion instability: Flow undergoes sudden, large 
amplitude excursion to a new, stable operating condition, due to the internal pressure gradient 
being less than the external pressure gradient (imposed by the pump). 

Geysering; an erratic boiling region exists, where the wall temperature moves between boiling 
and natural circulation in irregular cycles, due to the presence of gas; effect disappears at 
higher heat fluxes and higher pressures. 

Dynamic instability: a flow is subjected to dynamic instability when inertia and other 
feedback effects have an essential part in the process. Examples of dynamic instabilities are 
listed below: 

Density wave oscillations: Occur due to delay and feedback effects between flow rate, density 
andAp . 

Pressure drop oscillations: These oscillations occur due to dynamic interaction between 
channel and the compressible volume, which is initiated by flow excursions. 

Acoustic oscillations: These are high frequency (10-100 Hz) oscillations related to time 
required for pressure wave propagation; occurs due to resonance of pressure waves. 

Thermal oscillations: These oscillations occur in film boiling due to interaction of variable 
heat transfer coefficients with flow dynamics. 
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BWR instability: It is strong only for small fuel time constant and under low pressures; occurs 
due to interaction of void reactivity coupling with flow dynamics and heat transfer. 

Parallel channel instability: there are various modes of flow redistribution; occurs due to 
interaction between small numbers of parallel channels. 

The present work analyses geysering instability. 


1.4 Geysering 

1.4.1 A Brief Description 

The concept of geysering has evolved from the natural geyser. The mechanism of geysering 
was first described by Bunsen [1958]. He described the mechanism as follows: 

“If heat is added to a column of liquid which opens into a chamber of larger diameter, 
then it will ultimately rise to slightly above it’s saturation temperature. When this occurs, 
bubbles will form. Due to this bubble formation remaining liquid in the column will 
experience a decrease in hydrostatic pressure at a virtually constant temperature. Consequently 
more bubbles will form, further reducing the pressure, and ultimately much of the liquid in the 
column will be blown out. This phenomenon of blowing out is known as geysering”. 

Griffith [1962] suggested that two mechanisms, rather than the one given by Bunsen, 
could lead to geysering. The two mechanisms proposed by him are as follows: 

(1) Some amount of superheat in the liquid is required to nucleate a bubble. The superheat is 
often very large so that when a bubble forms, it grows very rapidly and often ejects most of the 
liquid which lies above it. In a nut shell, rapid growth of the bubbles is mainly responsible for 
geysering. 

(2) The second theory says that, even when the liquid is in a local saturation state then also 
some amount of bubbles will form and eventually lead to geysering. But disturbance of any 
size will not lead to geysering. For example, if the bubble formed grows slowly enough (small 
disturbance) such that it rises out of the column before any other can get started, then 
geysering may not occur at all. In other words, the decrease in pressure due to the formation 
of bubble should be large enough to form another bubble to start geysering. 
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Earlier researchers mostly concentrated on superheated liquid and saturated liquid for 
geysering. But early nineties experiments show that geysering also occurs at low pressures and 
low flows condition. The concern for possibility of geysering or condensation induced 
instability during start up from low pressure and low flow conditions in a natural circulation 
plant like SBWR (Simplified Boiling Water Reactor) was first raised by Aritomi et a/., [1992]. 
Then researchers shifted their attention towards sub-cooled liquid. The mechanism of 
geysering has been modified as follows: 

“Even though the bulk liquid is sub-cooled, the bubble at the vicinity of the fuel rod is 
superheated. Hence void will generate at some nucleating sites. At low flow rates and low 
pressure, the vapor bubbles produced in the heated channel will coalesce into a larger bubble. 
This will further grow due to decrease in hydrostatic pressure. This bubble will move towards 
the exit as its density is lower than the liquid. As the bulk liquid is still sub-cooled, the bubble 
will come in contact with the sub-cooled liquid and condense rapidly. Then due to collapse of 
bubble, pressure will drop suddenly. In order to balance this pressure drop, the sub-cooled 
liquid from the top will flow into the channel. Thus a flow reversal is induced. This process 
repeats periodically causing flow oscillation which is known as geysering”. [Aritomi et al., 
1992] 

1.4.2 Its Importance 

Natural circulation systems require power to initiate the circulation through void generation. 
This means the natural circulation reactor would be heated by fission energy from the startup 
under low temperature and low pressure condition. Thermal hydraulic instabilities have been 
reported under low pressure conditions [Chiang et al., 1994]. If thermal hydraulic instabilities 
were to occur at startup then the reactor would not continue operation during increase of power 
because the void fraction fluctuation in the reactor core would oscillate the reactivity. 
Therefore it is necessary to investigate and understand properly the role of reactivity during 
start up. 
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1.4.3 Interaction with Reactivity 


Reactivity is defined as the relative departure of neutron production factor (K^) from unity, 

, . _ . . Neutron Production Rate - Neutron Removal Rate 

that IS, Reactivity = 

Neutron Production Rate 


Mathematically, p = 



Neutron is required for the fission of nuclear fuel (sayt/^^). From every fission reaction, at an 
average 2.5 neutrons are evolved. If all the evolved neutrons react then the fission chain 
reaction will become uncontrollable. So we should allow only one neutron to react. That is we 
should make the reactor critical. But in actual practice the reactor is made slightly supercritical 
i.e. p = 1.01 to ensure that reactor remains critical, since in the reaction process some neutrons 
usually get lost. On the other hand if p < 1 then the number of neutrons available for the 
reaction will not be sufficient to move forward the chain reaction. Hence the reaction will die 
down and ultimately the reactor will stop. 

Hence, in a way reactivity can be used to control the reactor. As reactivity depends 
upon temperature and there is a mutual dependence between temperature and void fraction, a 
slight variation in void fraction would oscillate the reactivity. This makes the reactor difficult 
to control. 


1.4.4 Further Advances 

Aritomi et al. [1992, 1993] and Chiang et al. [1992, 1994] have conducted extensive research 
in the area of geysering under natural circulation. Inada et al. [1992] and Masuhara et al. 
[1993] have also performed small scale experiments to demonstrate this phenomenon. These 
experiments illustrated that the geysering mode oscillation would occur at low pressures and 
low flow conditions. 

Earlier attempts to model startup instabilities [Aritomi et al., 19^2] and [Paniagua et 
al., 1996] indicate that, to predict the possible startup instabilities correctly, it is important to 
accurately predict the vapor generation rate. In the present code a four equation drift flux 
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model has been used which is numerically more stable than the five equation and six equation 
model. Moreover, most of the models are envisaged for system with high pressures (greater 
than 20 bars) and thus are unsuitable for the simulation of geysering owing to the great 
influence of pressure on the void content in the sub-cooled boiling regime. Finally in the 
existing models there is no conclusive and physically well-defined description of mass transfer 
rate between the vapor and liquid phase, an aspect which is of great significance for 
understanding the sub-cooled boiling regime. In the present code the effect of the bubble 
formation as well as the condensation rate has been considered. 


1.5 Objective of the present work 


To study the two-phase flow phenomena, one can use various models such as Homogeneous 
Equilibrium Model (HEM), Drift flux model or a six-equation (two-fluid) model. HEM model 
assumes both the phases are in thermodynamic and mechanical equilibrium, i.e. they have 
same velocity and temperature. But geysering is a highly non-equilibrium phenomena because 
it considers growth of bubble due to sub-cooled boiling and its condensation in sub-cooled 
water in riser. So a HEM model cannot address this problem. Also in homogeneous model 
relative velocity between vapor and liquid is not considered. However, in two-phase flow there 
is always a relative motion between the liquid and vapor phase. 

The two-fluid models are the most robust and accurate. But their accuracies depend on 
the interfacial relationships for mass, energy and momentum transfer. The codes developed for 
these phenomenon, use empirical relationships which are not well established and they are 
continuously being upgraded with time. The second problem with two-fluid models is the use 
of single pressure for both the phases. This means that the velocity of sound in both the 
medium (liquid and vapor) will be equal. This makes the equations mathematically ill-posed. 
To avoid this problem, many codes use a numerical scheme to suppress the numerical 
instability. However, by this way, the codes miss the physical instability because the numerical 
scheme so used tries to damp the physical oscillation also. 
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The remaining possibility is the drift flux model. Its accuracy is better than a 
Homogeneous Equilibrium Model but lower than the two-fluid models. Because it considers 
the dynamic phasic slip between the phases by taking into account the vapor drift velocity iy ^ ) 

in the two-phase region. Earlier some codes based on drift flux model have been developed 
[RETRAN, ATHLET]. But these codes are applicable to forced circulation loop only. Also 
ATHLET code uses a five equation system that uses separate mass and energy conservation 
equations for steam and water, and a single mixture momentum conservation 
equation. [A.Hainoun et a/., 1996] This makes the system physically ill-posed as it has 
considered separate masses and velocities for steam and water, but still a single momentum for 
both the phases. 

In the present work a four equation drift flux model has been considered. It is more 
stable than the five and six equation models as described before. Most of the times one uses a 
semi implicit-explicit scheme which is better than the explicit scheme so far as time step is 
concerned and implicit scheme so far as numerical oscillations are concerned. Further the 
equations are mathematically well posed. The objective of this thesis is to develop and validate 
a code based on four equation drift flux model incorporating natural circulation loop features 
and with the proper consideration of vapor generation and vapor condensation model. 


1.6 Organisation of the present work 


The thesis is organized as follows; 

Chapter 1 gives a general introduction about the geysering phenomena and the literature 
survey. 

Chapter 2 describes the basis of Drift Flux model formulation with discretization equations 
and the Boussinesq approximation for the natural circulation loop. 

Chapter 3 deals with the simulation of small break LOCA using the two-phase code. 

Chapter 4 presents void generation and condensation model along with the McMaster 
experiment validation results. 
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Chapter 5 deals with the validation of the code with the Apsara Loop experimental data 
obtained from Bhabha Atomic Research Centre (BARC), Mumbai. 



Figure 1.1: Schematic of Natural Circulation Loop. 


14 




The Drift Flux Model 


Chapter 2 


The Drift Flux Model 

2.1 Drift Flux Formulation 
2.1.1 Governing Equations 


The Drift Flux model has been widely used for vertical flows in pipes and allows for slip 
between the phases and for a non-homogeneous distribution of the phases in a cross section. 
The basic model is due to Zuber and Findlay [1965] and contains elements of models by 
Wallis and by Bankoff [I960]. Its chief purpose is to provide a mechanistic framework for the 
calculation of the vapor void fraction . 

The drift-flux model for a two-phase mixture consists of two mass conservation 
equations, one momentum conservation equation, and one internal energy equation [Ishii, 
1975]. Here the following form of these four partial differential equations has been used: 


Total mass balance equation: 

dp„ ^ 1 ^ 

dt A dy 


Vapor mass balance equation: 


ap. 1 a 


■+ 


dt A dy 


p V + 


PgPt 


V, 


= r 


( 2 . 1 ) 


( 2 . 2 ) 
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Mixture momentum balance equation: 


dt A dy 


\ v^+BiElv^^ 

m ^ r 


dp . 

= + PmSy +/vi 


(2.3) 


Mixture energy balance equation 

PjnY^ + 

~ m m ni 




dt A dy 





y ^ PgPl 

r 1 1 1 

V 

A By 

Pn, 

[p; P:\ 

r 


+ KVr'+Kis+Q 


(2.4) 


The above four equations have been written using seventeen variables. The gravitational term 
g^, is a function of independent variable y only and does not depend on the state variables. 


Hence assuming that g^. is known from the geometrical configuration of the particular system 

under study, we must therefore specify twelve additional relationships among these variables. 
Two relationships may be obtained from basic definitions of mixture quantities. These 
definitions are; 

1. Definition of : 


p,„ =ap/+(l-a)p,‘* =p^ +p, 

2. Definition of /„ : 

m 


(2.5) 


( 2 . 6 ) 


2.1.2 Constitutive Equations 

The ten additional relationships needed for the description of two phase system are referred to 
as constitutive equations. A general form of each of these equations is presented below: 

1 . Thermal equation of state for the liquid: 

p,=p,{p,l,) (2.7) 
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This relationship must apply for a superheated as well as a sub-cooled liquid. 

2. Thermal equation of state for the vapor: 

( 2 - 8 ) 

This relationship must apply for both superheated and sub-cooled vapor. 

3. Caloric equation of state for the liquid: 

(2.9) 

This has been assumed that this relationship can be inverted to obtain 

I,=I,(p,T,) (2.10) 

4. Caloric equation of state for the vapor: 

( 2 . 11 ) 

We again assume inversion is possible, yielding 

( 2 . 12 ) 

5. Wall heat source: The heat source term in the energy equation is assumed to be of the form 

Q = hA^.if',-Tj) (2.13) 


where is the wall area, T^. is the wall temperature, and Tj- is the fluid temperature, 

possibly equals to T, or or some average temperature. The heat transfer coefficient h is 

obtained from correlations involving the fluid properties and the wall temperatures. 

6. Wall friction: 


fvis = q 


2d, 


(2.14) 


where c?,, is the hydraulic diameter and c is a friction multiplier depending strongly upon the 


void fraction a and other fluid properties. 

7. Relative velocity correlation: We assume the existence of a relationship specifying the 
relative velocity in terms of other state variables. This relationship will usually be strongly 
dependent on an assumed flow topology. 

8. Rate of phase change: The rate of phase change must be specified as a function of other 

variables. Various models employing equilibrium or non-equilibrium assumptions are 
possible. 
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9. Saturation curve: A saturation curve of the form 

T.=tXp) (2.15) 

must be specified. The saturation temperature provided by this curve enters the phase change 
expression above. 

10. A thermal constraint: The drift-flux model equations provide only a single mixture energy 
equation for the two specific internal energies 7^ andlj. Hence, an additional thermal 

constraint is necessary to partition the mixture energy into the liquid and vapor phases. A 
variety of such constraints is possible, the simplest of which is thermal equilibrium between 
phases. Another possibility is the assumption that either the vapor or the liquid is at saturation, 
depending upon whether vaporization or condensation is occurring. 


2.2 Drift Velocity 


The vapor drift velocity (v^ ) is the velocity of the vapor phase with respect to the velocity of 

center of volume or volumetric flux density of mixture. 

Mathematically, 

v,=v,-j (2.16) 


where j is the volumetric flux density of mixture. 

The mean volumetric flux density, or the superficial velocity, of each phase is defined as the 
volumetric flow rate of that phase divided by the total cross-sectional flow area. 

2,=—^ Bute, = A'',. = ( 2 . 17 ) 

Aj + A, + Ag 


Similarly j =■ 


2 . 




A + A + \ 

The mixture volumetric flux density will be given by: 

J = ji+J, ={^-ocy,+ccV^ 

Now, =v^-j=v^-il-ay,+ | 


(2.18) 

(2.19) 


18 



The Drift Flux Mode! 


= (l-aX''.-V,) = {l-aK 

where is the relative velocity of vapor with respect to the liquid phase. 


( 2 . 20 ) 


2.3 The Finite Difference Scheme 


The convective terms in equations (2.1)-(2.4) have been' written in a divergence or 
conservation form. Because there are large sources and sinks of momentum in a nuclear 
reactor (pumps, orifice etc.), it has not been attempted to ensure that the difference scheme be 
rigorously conservative in the treatment of momentum convection. Therefore, equation (2.3) 
has been rewritten in a nonconservative form that is more convenient for our applications. 
Multiplying equation (2.1) by and subtracting from equation (2.3) yields the equation 


dVn, X, dV„ 
2- + 1/ 2 

dt 


1 a 


dy P„A dy 


^PbPi yl 


1 dp 

= -T+gy 

Pm ^ 


Pm 


( 2 . 21 ) 


The advantage of using equation (2.21) over equation (2.3) is that rather than(pV")"^’ , is 
calculated directly. Its disadvantage is that it is not in conservation form, so we do not get 
rigorous conservation of momentum in the difference approximation. However, when the 
percentage changes in dependent variables from one cell to the next cell are not large, ±e 
nonconservation form of the momentum equation should not cause any problems. Equations 
(2.1), (2.2), (2.4) and (2.21) have been used throughout the following analysis. 


j-1 j-1/2 j j+1/2 j+1 

Figure 2.1: Mesh cell labeling convention. 
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V 


j-1/2 





p.p 


^j+1/2 


j-1/2 j j+1/2 

Figure 2.2: Mesh variables inside a computational cell. 


The mesh cell configuration and the labeling conventions for cell edges and cell centers are 
depicted in the above figures. The mass and energy equations are differenced over the mesh 
cells indicated by solid lines in Figure 2.1; the momentum equation is differenced over the 
dashed mesh cells. This forms a staggered spatial difference scheme. For finding the 
relationships among the variables at the edges and center of mesh cells, a “weighted donor 
cell” difference scheme has been used. The resulting difference equations are given below. 
From equation (2.1) the mixture density is updated as: 


ip.)T=ipJ, 


At 


AjAyj 


j*-. 


ip V ) . \ 

\r' m m / j-k— 



where 


( 2 . 22 ) 


(p.v, ),.i = \ k )j ip . ), + (p. U, }+ P\{V. )j|ip. )j - (p. 3 (2.23) 

(p.v. ),-l = 2 k )j., ip. ),_, + (p, ), }+ /3|(7. V, lip. - (p. ) J (2.24) 

Quantities on the right hand side of equation (2.22) are evaluated using n level values for p„ 


and the available n + 1 level values of V., . 

m 

Vapor density is updated from equation (2.2) as given below: 


(p, j*' =(p.);, +A,j-^[A^(p/,^-A_,(p.y.^]4(r.), 


(2.25) 
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where 


ifi.y. u = I k I - (p. I K 4''. Ih I - 1 3 


(2.26) 




(2.27) 


Quantities on the right hand side of equation (2.25) are updated using n level values of (p^ 
and the available n + 1 level values of (v^ . 

The finite difference equation used to approximate equation (2.21) is 


(v.r\ -(yJ. . 

At 


+(yj’ 




\yX,-(yJ, 


Ay 


{PmT. 1^. 





n 

'Ap^pyy 

n 

L 

y+1 


j ^ 


(pj. 


^ pTXpT' ^ 

Ay 




+ (^r),>i + 
2 


«/ VIS 

\ jj*- 


(2.28) 


- CFL UX - DFLUX + ^ + g,+ 

p, ,(Ay^+Ay^,,) 


(^J;=(^J,+Ar 

)j is the explicit estimate of (V„ . 


^^2 


^ VL 

V JJ 


(2.29) 


where 


(p ). , Ay^.^p.+Ay^.p^.,, 


2 Ayj +Ay 


2+1 


(2.30) 
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CFLUX = K 


ay 1 

m 

dy ~ 1 


+ 


<y.)j 


(V. )j - (V'. ^ (t'. V, - fy. ), 




4Vj>, 






^yi 




(2.31) 


DFLUX 


1 a 


Pm^dy 


^P.Pl y2 


p.^Mj + ^j^il^yj+^yjJ 






\ j j+i 


{(Vjj 


+ 


+ 


^ |(V. ), + ),., |[(V, ), - (V, ),., I- A / 1(V, V, + (V, I } 


Pn, 


Jj 




(2.32) 


Finally, 


2Ar 

P. iWy +Ay;^,) 


[Pj^x - 


P%i-Pj + P 


;) 


'"2 


(2.33) 


The energy equation (2.4) is finite differenced as: 


(n I (n T Y ^ \{PmK^m)i+l~ ^ . liPm^m^m)j-l 


At 


^^yj 


+ 


1 

AjAyj 



n 

'PlP,{h-h)y' 

^ r 

n 

L P". 

. \ 

^^2 


. 1 
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+ [kv;1+[wJ]+q". (2.34) 
Mixture energy equation (2.34) is rearranged as 


(/Jr +At{-FLUXG-FLUXL-PWK 


+ 'f; ((K ); +(V',V, )+(«', *),+Q;} 

(2.35) 

where 


“ A L V u - u 

^Ayj L ^ 2 2 ^2 2j 

(2.36) 

FLUXL^ \ A ,(p,;,V,) , A ,(p,/,V,) _, 

(2.37) 

The convective energy flux terms are given by 



(2.38) 

2 

(2.39) 

{p,I,V,)j,l =|feU(p,/,)j +(P,/,)j..]+^|V',|[(P,/,), -(P, /,);.,]} 

2 jiL 

(2.40) 

(P,I,V,)jJ. "ItKViiP,/,);-, +(p, /,),]+ ^M[(P,y,)j-, -(P,',)J 

(2.41) 
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The vapor and liquid velocities used in these expressions are defined as 




1 



Ilk. 

(Pm^i 

-) 


;+ 



(2.42) 

(2.43) 


The pressure work term is approximated as 

’(p*Li-(pJ;4+P° (P*U 
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f r 
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PWV — ^ 
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^ 2 2 
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(VrU 


(2.44) 


Quantities on the right hand side of the energy equation (2.35) are evaluated using n level 
values for p„,Pg and/^, whereas n + 1 level values are used for (y„)j,Pj and the overall 

divisor These difference equations have been supplemented by additional relationships 

among variables at the edges and center of mesh cells. The “weighted donor cell” difference 
scheme that has been used above is given by: 


a , = 


LLi. 

K 2 , 


pj + 




Pi- 


j+i 


\OCj + 


1-/3 
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^ 1 =1 
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fi-P] 
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[ 2 J 
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1 2 J 




(2.45) 

(2.46) 

(2.47) 
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(2.48) 


where /3 is the weighting parameter. It varies between -1 and +1. It can be chosen in various 
ways. For example, 

jS = 0 yields more accurate but less stable central difference technique. 


/3 



V , = ±1 gives full donor cell difference or fully upstream or fully upwind 

/ 


difference scheme. This scheme is stable provided fluid does not convect through more than 
one mesh cell in one time step. Though this technique is stable, it is first order accurate. 


2.4 The Solution Procedure 


The numerical scheme that has been developed is based on SIMPLE (Semi-Implicit Method 
for Pressure Linked Equations) algorithm formulated by Patankar [1980]. Primitive variables 
p, V, and I (ultimately it will lead to temperature T) are used in a staggered grid system. The 
computation domain is divided into rectangular control volumes with one grid point located at 
the center of the control volume that forms the basic cell. Pressure p and temperature T are 
calculated at the center of the cell. Velocity V is calculated for points that lie on the face of 
these basic cells. A typical basic cell with mesh cell labeling convention is shown in the 
following figure. 



j-1 j-1/2 j j+1/2 j+1 
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The basic philosophy of the pressure correction method used in SIMPLE algorithm is as 
follows; 

1. Start the iterative process by guessing the pressure field. Denote the guessed pressure by 
Po • 

2. Use the value of to solve for the velocity V in the momentum equation. Since this 
velocity is associated with the values of p^, denote this by Vq • 

3. Since the velocity is obtained from guessed value of pressure the velocity V'g, when 
substituted into the continuity equation, will not necessarily satisfy that equation. Hence, 
using the continuity equation, construct a pressure correction p' which when added to Pq 

will bring the velocity field more into agreement with the continuity equation. That is, the 
“corrected” pressure p is 

P = Po+P' 

Corresponding velocity correction V' can be obtained from p' such that 

V=Vo+^' 

4. Designate the new value p as the new value of pg . Return to STEP-2, and repeat the 

process until a velocity field is found that does satisfy the continuity equation. When this is 
achieved, the correct flow field is at hand. 


2.5 Stability 


One has to use a proper time step and grid size to avoid the problem of numerical difrusion. 
Numerical diffusion means as we discretise a differential equation by FDM (finite difference 
method), normally we neglect the higher order terms. Suppose they diffuse from one 
calculation to the other and the error keeps on building, then it may cause numerical error. As 
we know from Von-Neumann stability criteria that explicit method is conditionally stable and 
has got a stringent time step limitation. Using a smaller time step means number of iterations 
to reach at the solution is more. This adds to the round off eiror. On the other hand, implicit 
methods are unconditionally stable. But it uses the advanced time step value to proceed to the 
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next cycles, which are generally not known. Hence a semi implicit-explicit scheme has been 
used which is better than the explicit scheme as far as time step is concerned and implicit 
scheme so far numerical oscillation are concerned. Basically the iterative explicit treatment of 
pressure, equation of state term makes this technique applicable at all flow speeds (sonic as 
well as sub sonic). 

The necessary condition for stability is that the fluid should move less than one mesh 
cell per iteration. This is known as effective Courant number criteria. Hence the stability of the 
numerical scheme can be assured only for the time step sizes below the effective Courant 
number as given by: 


Ay 


<1 


Also, stability is expected when the upwinding weighting parameter /3 is chosen such that 


P 


> maxi 


Ay 


pJ^y 


In the present work effective Courant number has been taken less than 0.25. Also the grid size 
should be chosen so that it is larger than the bubble radius. Otherwise there may be instability 
due do sudden change in state in two adjacent cells due to the presence of bubble in one cell 
and liquid in the next cell. On the other hand, too large a grid size will make the truncation 
error prominent. 


2.6 Closure 

In this chapter, formulation of Drift flux equation is described along with the finite difference 
discretized form for these governing equations. The solution procedure for solving the system 
of equation has been described. This chapter also presents the general concept of stability 
criterion. 
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Chapter 3 

Simulation of Small LOCA 


3.1 Motivation: The Three Mile Island accident 


The motivation behind choosing this problem is to illustrate how a small break LOCA can 
trigger to a large scale accident. In March 1979, an event occurred at the Three Mile Island 
Unit 2 (TMI-2) that resulted in the first case of melted fuel in a full scale commercial nuclear 
reactor power plant. The initial phase of this accident involved a loss of steam generator feed 
water that led to depletion (boil-off) of the liquid inventory on the steam generator secondary 
side. The resulting loss of heat removal caused heat-up and expansion of the primary coolant. 
The pressure rose in the reactor primary cooling system until the reactor shutdown. This led to 
the actuation of the pressurizer relief valve. This valve failed to close when it was supposed to- 
after pressure dropped below the set point for closure. This initial stage of the incident does 
not call for sophisticated analysis, and the models employed in the plant simulators for 
oj)erator training could have addressed it adequately. 

That, however, was not the case due to the occurrence of the subsequent small break 
LOCA. The simulators available at that time did not include provisions to account for the 
presence of two-phase flow within the primary coolant system. Any void (steam) generated 
anywhere in the system was thought to end up, without delay, in the upper portion of the 
pressurizer. The simulators therefore taught the operators to concentrate on keeping the 
pressurizer liquid level within prescribed limit; a liquid-full pressurizer was regarded as an 
indication of a liquid-full primary coolant system. 
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Lack of heat removal combined with loss of pressure caused by discharge of coolant 
through the stuck-open relief valve, eventually caused liquid flashing in parts of the reactor 
vessel and accumulation of steam within its upper regions. The liquid displaced by the steam 
was pushed into the pressurizer, giving the high water level indication and causing operator to 
turn off the emergency water injection pump. This ultimately leads to excessive fuel heating 
and fuel melt down 


3.2 Problem Statement 


A schematic of the problem geometry is shown in Figure 3.1 A constant pressure pump 
( p =75 bar) induces flow of a liquid (water) at a temperature of 500°K through a single 
branching system into a pressurized vessel at a pressure of 70 bar. The flow is calculated until 
steady state is achieved. A guillotine break is then made in one of the line branches at its 
vessel connection and the transient output is then calculated. Figure 3.2 shows the dimensions 
and computational zone for the pipe system, which contains three components and one 
junction. The bottom end of the component (pipe) 1, which is connected to the vessel, is 
treated as an isolated end with a constant pressure outflow boundary condition. At the time of 
the break the boundary condition of bottom end of pipe 1 will be changed to ambient condition 
{ p =lbar) to reflect the pipe’s separation from the vessel. The top end of the pipe 1 and the 
bottom end of the pipe 2 and pipe3 are connected to the junction. The top end of the pipe 2 
remains connected to the vessel during the entire calculation tind is described by an isolated 
constant pressure outflow boundary. The top end of pipe 3 is treated as an isolated constant 
pressure inflow boundary to simulate the pump. 
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3.2.1 Geometry 



Figure 3.1: Single loop system geometry for the test problem. 


3.2.2 Computational Domain 


OUTLET 


25 H 
22 

20 


300 cm 


A 

L 

— 5y=5.25 cm 


_5y=22.5 cm 


5y=23.75 cm 


c 


15 cm dia 


30 








38 


INLET 


200 cm 


7 

5-t 


- 10 cm dia 


- 7 cm dia 


A: Connected to Vessel before and after break 
B: Location where the break occured 
C: Connected to Pump side 


OUTLET 

Figure 3.2 : Structure of the computing mesh for the test problem. 
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3.3 Constitutive Relations 

Basic definitions of mixture quantities are as follows: 


p,={\-a)p° 


P»=P,+P/=«p/ + (l-«)p/‘’ 


a = 


Pi -Pn,+P, 

p/° 


1. Thermal equation of state for the liquid: 

p = Po + a'(p-p,°) fora<a^ 


2. Thermal equation of state for the vapor: 

p = (y-l)p fora > a. 


3. Caloric equation of state for the liquid: 

I,=E,+C,{T,-T„)-C„{T,-Tj 

4. Caloric equation of state for the vapor: 


5. Wall heat source: 

Q = Heat flux density = 


Power 


rzRjX Length of heated section 


6. Wall friction: 

The distributed losses (pipe wall friction and local losses that occur 
sudden change in area) can be given as: 


(3.1) 

(3.2) 

(3.3) 

(3.4) 

(3.5) 

(3.6) 

(3.7) 

(3.8) 

(3.9) 

due to 
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/. 

J VIS 


=_x 

R. 


-|2 


Pi 




pyj-jpyj 


Using Govier and Aziz [1972] correlation, the friction coefficient / is given by: 


f = a + 


Re" 


where, 


a = 0.026{k/2 R, +0.1 33(fe/2i?, ) 
b = 22.0{k/2Rj-^ 
c = l.62{k/2Rj' 


0.134 


R^ = Hydraulic radius 


(3.10) 


(3.11a) 

(3.11b) 
(3.11c) 
(3.1 Id) 


(t)j.p = Two-phase friction multiplier 

The coefficient relates to the local losses and is given by 
f(l I'lR 1 

^ JAJ — U. Where {L/2R^ ) is the hydraulic diameter of an equivalent straight channel 

Ay 

and Ay is the segment length over which the loss occurs. 


7. Relative velocity correlation: 


y/ accounts for the relative velocity effects and is given by: 


¥=P, 


1 + (P«-Pj 


P V" 

r^m m 


(3.12) 


8. Rate of phase change: 

The vapor generation rate is calculated as follows; 

i 

q-S 




q = Interfacial heat flux 


evap 


h^^ = Latent heat of vaporization 


(3.13) 
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^ I 

Hence, 

r 


(3.14) 


(3.15) 


Where, Z= Thickness of thermal boundary layer over which the liquid temperature 
changes from its interior bulk temperature T, to the value 

Theofanous et al. [1975] had given the value of I = l„, for a single, nontranslating bubble 
growing in an infinite fluid region as; 


6 p; C,\T,-T, 


(3.16) 


Moalem and Sideman [1973] had given the value of 1 = 1, for bubbles which are translating 
with respect to surrounding liquid with speed V , as 




K 


V/2 


V j 


(3.17) 


Where, RCj, = 


2rVp° 


= Bubble Reynolds number 


Pr = 



= Liquid Prandtl number 


(3.18) 

(3.19) 


Combining both of these effects. 


1 

/ 



L I, 


(3.20) 


9. Saturation curve: 

r, =255.2 + 117.8/7 


(3.21) 
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10. Thermal constraint: 

For equal-phase temperatures, the mixture temperature can be computed from the 
mixture internal energy as the solution of a quadratic equation 

Pm^m ~ Pg^g PJi (3.22) 

When the vapor is considered to be saturated, its temperature is determined from the 
mixture pressure by the relation 

=255.2 + 117.8p‘'-^-^ (3.23) 


The Momentum exchange function: 


K = 


pS 


8a, 




I2v 


Where, 


a, = a , V = V; (l - a)'^’^ fora < 0.5 

a, =l-a, v=v^a~''^ fora >0.5 

The surface area per unit volume of bubbles with mean radius q, is given by: 


(3.24) 


(3.25) 

(3.26) 


S = — fora <0.5 
ro 

S = fora >0.5 


(3.27) 

(3.28) 


The mean bubble radius is given by: 
( 3a 


^0 = 


4;iA^ 

/ 

' 3(1 -g) ' 
AnN 


1 / 3 


fora < 0.5 


fora > 0.5 


(3.29) 

(3.30) 


The maximum stable bubble radius in terms of critical Weber number and o , the interfacial 
surface tension is. 



(3.31) 

’2pX' 

Mm(ro,r„,) 

(3.32) 
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Relative speed used in equation (3.31) is given by, 



= Mass averaged mixture velocity = -^j [(V ,„ ). + (V,„ ] 

= Average relative velocity = [(V ^ ). + (y ^ ] 

P, = Parameter accounting for turbulent fluctuation = 0.1 


(3.33) 


3.4 Mesh Construction 

3.4.1 Mesh construction for the components 





JT3 

JT2 

JTl 


JBO 

JBl 

JB2 


j 



Figure 33: Arrangement of fictitious cells at the bottom and top of the component ‘I’. 

Referring to Figure 3.3 
I: denotes component (pipe) 

J: denotes cell 

JBO: denotes first fictitious bottom cell 

JB 1 : denotes second fictitious bottom cell 

JB2: denotes first real bottom cell 

JTl: denotes last fictitious top cell 

JT2: denotes last real top cell 

Component can be divided into different segments. 

Assumption: All cells in a segment has identical properties. 
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If necessary, we can define different properties for every cell by making each cell a segment. 
The first two bottom cells JBO and JB 1 and last top cell JTl in the component are “dummy” or 
“fictitious cells used to set boundary conditions and accomplish coupling of components. 
Two cells are needed at the beginning of a component to set the boundary conditions because a 
staggered mesh arrangement has been used, where velocities are located at the cell boundaries. 
All other dependent variables are located at the cell centers. 


3.4.2 Mesh Construction for a Junction Cell 

A junction cell is a three-dimensional rectangular control volume that can join up to four 
components in a plane. In the present case flow in the third dimension is restricted. This 
restriction allows the use of a two-dimensional rather than a three-dimensional solution 
algorithm for the junction cell. 

A one-dimensional component may be attached perpendicular to any of the four faces 
of the junction as shown in the Figure 3.4 and Figure 3.5. For simplicity only one component 
is allowed per face and components opposite to one another must have the same cross- 
sectional area at the junction. If no components are attached to a given junction face, it is 
treated as a rigid wall boundary. 

The coupling of a junction cell and two or more component is accomplished by making 
the junction cell the neighboring cell to the JB2 (first bottom real cell) or JT2 (last top real 
cell) cells that belong to the adjoining ends of the components. The JBl or JTl fictitious cells 
at the adjoining ends overlap the junction cell. The centers of the overlapping cells coincide 
with the center of the junction cell so that cell centered quantities in the junction cell and cells 
JBl and JTl are identical. Furthermore, the boundaries of the overlapping cells on which 
velocities are stored coincide with the boundaries of the junction cell so that velocities on the 
common boundaries are identical. 
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LANG(I)=1 

KTOP(I)=K 


LANG(I)=4, KTOP(I)=K 


JT2 


JT2 



JT2 


JT2 


LANG(I)=2, KTOP(I)=K 


LANG(I)=3 

KTOP(I)=K 


Figure 3.4: Junction cell orientation for component coupling connecting from the top. 
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LANG(I)=1 
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Figure 3.5: Junction cell orientation for component coupling connecting from the bottom. 
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3.5 Component Boundary Conditions 


Various types of boundary conditions may be used at the ends of the one-dimensional 
component meshes. The following is the list of boundary conditions that are available in the 
present code. 

1. Prescribed Zero Velocity. A prescribed zero velocity can be set as boundary 
condition in the component. 

Example: V(750) = 0 

2. Prescribed Non-zero Velocity. A nonzero or time varying prescription can be set as 
boundary condition in the component. While setting nonzero velocity, one must also 
specify the pressure, void fraction and temperature at the boundary. 

Example; V{jB0)= constant or /(^) 

3. Uniform Outflow. A uniform or gradient free outflow boundary condition at the ends 
of the component can be specified by setting the velocity at a particular cell equals to 
the velocity at the previous cell. Use of this boundary condition also requires the 
additional specification of pressure, void fraction and temperature at the boundary. 
Example: V(75l)= V(/Sl-l) = V(y50) 

4. Periodic Coupling. The top and bottom boundaries of component can be joined 
together through the periodic boundary condition. 

Example: V(jB0)=V(jT2) 

5. Prescribed Pressure. When pressure is prescribed at the JBO end, it must be set for the 
cell JBl. The velocity V(JB0) should be set equal to v(jBl) so that fluid can flow 
freely into or out of the specified pressure region. The values for pressure, void fraction 
and temperature at the boundary also need to be specified. 

Initial and boundary conditions: The initial and boundary conditions used in the present 
work are given below. 

I.C.: at r = 0 , T= specified, p = specified, ct = 0 for all cells. 
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For f > 0 , B.C, 1: at A, B, C in Figure 3.2 T = specified 
B.C. 2: at A, B, C in Figure 3.2 p = specified 
B.C. 3: at A, B, C in Figure 3.2 cc = specified 

3.6 Solution Algorithm 

3.6.1 Solution Algorithm for Components 

A calculation cycle has been broken down into four tasks. First, the momentum equation is 
advanced explicitly using the previous cycle values. Next, iteration is made to replace the 
pressure used in the first task with advanced time values. Iteration is needed because the 
advanced pressures depend on the velocities being calculated. This part of the cycle contains 
the main implicitness of the numerical scheme. The third task in a cycle is to update all other 
dependent variables. Finally, the fourth task consists of data output and time step control. 

3.6.2 Solution Algorithm for a Junction Cell 

Junction cell values are calculated using the two-dimensional analogues of the one 
dimensional equations of motion [equation 3.1- equation 3.4]. The difference equations for 
junction cells are the two-dimensional counterparts of the one-dimensional finite difference 
equations. In the cross convection of momentum two-dimensionality is not maintained, 
because a one-dimensional component can have momentum only in the axial direction. This 
means, momentum entering the right or left side of a junction is not allowed to convect into 
components connected to the top or bottom of the junction. To ensure the coincidence between 
-the junction and fictitious cell centers and boundaries, an additional segment is added 
automatically to the end of the component that adjoins the junction cell. This additional 
segment contains the required one or two fictitious cells, depending on which end of the 
component it is added to, and one real cell. Inclusion of this real cell, however, slightly 
increases the overall length of the component. But a better approximation of the actual length 
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can be obtained by specifying input data that are one cell short of the desired length (two cells 
short if junctions exist at both ends of the components). 


3.7 Numerical Scheme 


The numerical scheme that has been developed is based on SIMPLE (Semi-Implicit Method 
for Pressure Linked Equations) algorithm formulated by Patankar [1980]. Primitive variables 
p, V, and I (ultimately it will lead to temperature T) are used in a staggered grid system. The 
computation domain is divided into rectangular control volumes with one grid point located at 
the center of the control volume that forms the basic cell. Pressure p and temperature T are 
calculated at the center of the cell. Velocity V is calculated for points that lie on the face of 
these basic cells. A typical basic cell with mesh cell labeling convention is shown in the 
following figure. 


j-1 j-1/2 j j+1/2 j+1 

The following is a brief description of the numerical procedure: 

1. Set or estimate the initial values of dependent variables: Vq.Pq and 

2. Advance in time {t + At) 

3. Evaluate the pseudo velocities from the momentum equation using the Vq.Pq and Tq 
value at n‘^ time. 

4. Use the pseudo velocities in continuity equation to find the updated density p andct . 
Also using this pseudo velocity, calculate the correction in continuity equation “D”. 
Our aim is to calculate such velocities for which “D” equals to zero. Actually D is an 
approximation to cell volume change per unit volume. This is given by, 
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D = 


A 9v 


AV.. = 


1 A 


^yj 


‘^,.1 - A . 


5. If “D is not zero, then using this D solve the pressure correction equation. 
F 


Ap = - 


'' 

V Jj 


where, = p-/(p„„p, J „)=0 . 

The equation of state /(p„ , p^ J,. ) is evaluated using 

p.=(pJV(i+o) 


and 


h=iiJ- 


(pj 


D 


6. Using this Ap solve the pressure field and velocity field as follows: 
p = Po+Ap 

Vj=Vj+ ^ 


V =v ~ 


2Ar * Ap 


p ,(Ay^._, +AyJ 


7. Return to step-3. Solve the momentum equation based on new pressure field p and 
continue updating values until the correction in continuity equation “D” satisfies a 
predetermined convergence criteria. 

The general form of convergence criteria is given by: 

|p-/fc.P«>^«)| 


< e ; Here e is taken to be 0.001. 


P + , I ^ )j 

8. Solve for the new temperature field using new velocity field values. 

9. When the convergence criterion for “D” is satisfied, we will have the value of primitive 

variable for the next time step i.e. V"*\p''*\T''*^ 

10. Return to step-2 and continue till steady state is reached. 
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3.8 Variable Time Steps and Subcycling 

One of the advanced features of the current code is the use of variable time steps and 
subcycling during calculations. This is due to the reason that network systems often contain 
low speed flow with slowly varying properties in one region and high speed flow or flow that 
requires a finely detailed description in another region. The variable time stepping and 
subcycling provisions in this code are designed specifically for such systems to provide 
efficient and accurate calculations. The time steps are determined by numerical stability 
requirements. 

For numerical stability the time step in each component and junction is limited by the 

Courant stability criterion, which is given by \ ~ • The time steps determined 

Ay{Ax) 4 

according to this criterion are then increased or decreased by 1%. The direction of adjustment 
is determined by the relative ease of the previous time integration of the system. If fewer than 
five system iterations were required, the time steps are increased; otherwise, they are 
decreased. The time step for the system is related to the maximum time step determined for the 
components and junctions. The system time step is not a computational time step, but rather 
provides a time level toward which all the components and junctions are integrated 
simultaneously. The system time step is increased or decreased by 1%, based on relative ease 
of the previous system integration and the number of subcycles used. The system time step is 
increased if the number of subcycles is less than a specified number of subcycles, otherwise it 
is decreased. 
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3.9 Results 



0 1000 2000 
Time in ms 


Figure 3.6: Time Variation of Pressure at different cells (break at 2000 ms). 
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Cell Number 


Figure 3.7: Time Variation of Pressure in the first pipe. 
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Figure 3.13: Time Variation of Velocity near the entry from the pump. 
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Figure 3.14: Enlarged View of Time Variation of Velocity near the entry from the pump. 



Figure 3.15: Velocity vector before the introduction of break 
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Figure 3.16: Velocity vector after the introduction of break. 
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Figure 3.17: Time Variation of Relative velocity at the break. 
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Figure 3.18: Time Variation of Relative velocity at the middle of the first pipe. 



Figure 3.19: Comparison of Relative velocity in the first pipe just after and long after the break. 
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i.22: Time Variation of Mass Flux at the top end of the third pipe. 


Mass Flux Bottom Pipe 1 

Mass Flux Top Pipe 2 



Time in ms 

Figure 3.23: Time Variation of Mass Flux in the whole system- 
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Figure 3.24: Time Variation of Density before and just after break calculated with single precision 
accuracy. 



Figure 3.25: Time Variation of Density before and just after break calculated with double precision 
accuracy. 
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Figure 3.27: Time Variation of Pressure, Mixture Density and Void Fraction before and long after break 
calculated with double precision accuracy. 
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Figure 3.29: Time Variation of Void fraction before and just after break calculated with single precision 
accuracy. 



Figure 3.30: Time Variation of Void fraction before and Just after break calculated with double precision 
accuracy. 
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Figure 331: Time Variation of Quality before and after the break calculated with single precision 
accuracy. 



Figure 3.32: Time Variation of Quality before and after the break calculated with double precision 
accuracy. 
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3.10 Discussions 


The variation of pressure at different cell of the component with respect to time has been 
plotted in the Figure 3.6. The flow is calculated up to 2000 ms where steady state is achieved. 
At this instant a guillotine break is made in one of the pipe at the vessel connection. The plot 
shows an instant fall in pressure at the first cell (Cell No. 2) of the broken pipe. After this point 
the pressure in the pipe fluctuate with time as expected. This is because the pressure at the cell 
suddenly reduces to a lower value (from 72 bar to 1 bar) and due to which flashing of steam 
starts. Due to production of bubble the pressure in the cell rises, which in turn reduces the 
number of bubble being produced. Hence again the pressure falls due to lack of bubble 
production and this process repeats. 

Figure 3.7 shows the variation of pressure in the first component at a particular time. In 
the present case pipe break is made at the lower vertical pipe. Hence the variation of pressure, 
before and after the introduction of break, is more in the vertical pipes than the horizontal pipe. 
Figure 3.8 and Figure 3.9 shows the variation of pressure in the second and in the third 
component respectively. The variations of mixture velocity with respect to time at /different 
cells are shown in Figure 3.10 - Figure 3.14. The square of sonic velocity in the liquid at this 

particular pressure (Ibar) is nearly 12340 That means the sonic velocity is 

around«l 1 \mls . Referring to the Figure 3.10 the velocity near the break is =10cm/m5 which 
equals to near sonic velocity. The fluctuation of velocity is due to the fluctuation of pressure 
governed by Bernoulli’s equation. Figure 3.15 and Figure 3.16 gives the vectorial 
representation of velocities in the component before and after the introduction of break. The 
variation of relative velocity after introducing break is given by Figure 3.17. The variation of 
relative velocity in the first pipe is shown in Figure 3.18 and Figure 3.19. This confirms that 
voids are created only after break is made. Figure 3.20, Figure 3.21 and Figure 3.22 gives the 
mass flux variation with respect to time in first, second and third pipe respectively (Refer 
Figure 3.2 for the pipe number notation). Figure 3.23 shows the mass flux variation in the 

whole system. The area under each curve in Figure 3.23 represents the mass flux in gmf cm} 
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from each of the component. The mass flux coming from pipe 3 goes to pipe 1 and pipe 2 and 
this can be verified in Figure 3.23. The effect of single precision and double precision 
accuracy calculation on the mixture density before and just after the break has been shown in 
Figure 3.24 and Figure 3.25. Here double precision calculation tends to give an erroneous 
result as the mixture density at one of the cell is even more than the single phase density. This 
may be due to numerical diffusion occurring in double precision calculation. The time 
variation of pressure, mixture density and mixture void calculated with single precision and 
double precision accuracy has been shown in Figure 3.26 and Figure 3.27. Figure 3.28 shows 
the time variation of vapor production rate after the introduction of break. The void fraction 
variation calculated with single precision accuracy and double precision accuracy has been 
shown in Figure 3.29 and Figure 3.30 respectively. Figure 3.29 and Figure 3.30 shows that the 
void fraction variation occurs only after the break and the fluctuation is more prominent in the 
cells of the broken pipe. The fluctuation of pressure changes the saturation temperature of 
liquid in that particular cell. This in turn varies the density of fluid in the cell and finally the 
density is responsible for varying the void fraction. Finally the variation of quality in the entire 
system at different time has been shown in the Figure 3.31 and Figure 3.32. 


3.11 Closure 


The simulation of small break LOCA on a test problem has been performed successfully and 
the various results that have been derived are quite encouraging. It was the motivating factor to 
include the void generation and condensation model in the present code, to be applied to sub- 
cooled boiling, which is described in the next chapter. 
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Chapter 4 


Application of Generation and Condensation Models 
in Sub-cooled Boiling 

4.1 McMaster Experimental Set-up 


This experimental project was jointly funded by NSERC (National Science and Engineering 
Research Council, Canada) and AECL (Atomic Energy of Canada Limited).The objective of 
this project was to obtain detailed measurements of void generation and collapse, void profiles 
of different bubble size and condensation rates during low-pressure sub-cooled boiling. The 
test section was divided into two parts: a heated section followed by an unheated section. 
Measurements were taken in both areas. In the heated section, the void results from the net 
effect of generation and condensation, whereas in the subsequent section, condensation has 
been studied as a separate effect. 

The test section was installed in a low pressure loop. A pump had been used to 
circulate the water to the loop through a preheater, which controls temperature to the test 
section. The test section was a vertical concentric annulus with measurements shown in the 
Figure 4.1. The annular channel consisted of an inner and an outer tube. Heat was generated 
uniformly in the test section by direct heating from a 55 kW stabilized direct current (dc) 
power supply. The heat was uniformly input via the inner tube (diameter 12.70 mm). The 
width of the annular channel was 6.15 mm. The test section had upward flow and consisted of 
a heated lower section 300 mm in length and a subsequent unheated upper section 120 mm in 
length. Flow was measured by a rotameter and void fraction by a traveling gamma 
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densitometer. In addition, high-speed photography was used to capture bubble action in the 
generation and condensation regions. The details of experimental procedure and void profile 
results, together with error estimation, were given by Donevski and Shoukri [1989]. The 
experimental parameters have been given in Table 4. 1 . 



Case 

1 

2 

3 

4 

5 

6 

Heat flux (W.m'^) 

714403 

586420 

594136 

576132 

481482 

732700 

Pressure (kPa) 

165.9 

211.4 

173.5 

168.7 

154.2 

152 

Mass flux (kg/m^s) 

367.5 

315.1 

420.0 

429.1 

392.1 

450 

Inlet temperature(°C) 

85.0 

98.1 

95.5 

95.2 

93.8 

85.8 

Inlet subcooling(°C) 

29.2 

23.8 

20.2 

19.6 

18.6 

26.3 


Table 4-1: Experimental parameters used in McMaster experiment 
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Figure 4.1: Geometry of the McMaster experimental test section. 
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4.2 Mathematical Formulation 


Liquid mass balance: 



Vapor mass balance: 


±r 

dt 


ap 


8 


+ 


dy 


ap^V 
8 8 




(4.2) 


Mixture momentum balance: 


dt 


(p V ) 

^ m m' 


a { 


dy 


8 8 


(l-a)pj 


Oy 2 




dp 

■^ + P 8 + / . 
dy tn y vis 


(4.3) 


Mixture energy balance: 



= KV'^ +W +0 
r vis 


(4.4) 


4.3 Modeling of Void Formation in the Sub-cooled Boiling 
Regime 


In the sub-cooled boiling regime, boiling occurs at the liquid cooled heating surfaces due to 
high heat flow densities, although the fluid has on average not yet reached the saturation 
temperature associated with the system pressure. In accordance with a suggestion by Griffith 
(Griffith et ai, 1958) four zones of heat transport and flow activity can be differentiated along 
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a channel axis until saturation boiling has been reached. Zone I is referred to as single phase 
heat transfer zone. Zone II is onset of nucleate boiling (ONB) zone, Zone HI is onset of 
significant Void (OSV) and Zone IV is the saturated core flow region. A visual of all these 
four zones is given in Figure 4.2. The equation of the following effect must be formulated in 
order to determine the void content in the sub-cooled boiling regime: (1) on set of nucleate 
boiling (ONB). (2) bubble formation and bubble growth (bubble generation rate) at the heating 
surface (3) onset of significant void (OSV) and bubble departure diameter (4) bubble 
condensation in the sub-cooled core flow. 



Figure 4.2: Temperature and Void fraction variation in sub-cooled and bulk boiling regions. 


64 




Application of Generation and Condensation Models in Sub-cooled Boiling 


4.3.1 Onset of Nucleate Boiling (ONB) 


It describes the point at which first boiling nuclei are activated at the boiling surface. A 
correlation by Bergels et al. [1981] and Rohsenow has been used to determine the wall 
superheating at which ONB will be activated. 


^^sup 


4 -1.156 

P 

1100 


0A63p 


0.0234 




r " 

4 -1.156 

P 

1100 


0.463p0-0234 


(4.5) 


(4.6) 


The range of validity of the correlation is: channel diameter 2.4 and 4.6mm; flow velocity: 3- 
7ms' ’.Pressure: 1-36 bar. 


4.3.2 Onset of Significant Void (OSV) 


A schematic description of heat transfer mechanism in sub-cooled boiling is given in Figure 
4.3. Assuming that the sub-cooled boiling occurs, the total heat Q supplied to the fluid from 
the heating surface per unit area can be arbitrarily partitioned into two fractions- single phase 
convection and heat required for growth of bubbles: 

Q. — (2i* 4 Q,gv 
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Heat required to generate void is composed of: Q^y = 

Xhc 3,t)Ovc two expressions 3.re combined, to ^ve the hc3.t flux 3s* 

It ti •• 

? “ ^evap (4 7 ) 

The heat flux associated with the single phase liquid convection from the heating surface is 
given as: 


Where T^. and 7] are the wall and liquid temperature respectively and /zj<j,is the heat transfer 

coefficient of the single phase liquid flow, which can be determined from the Dittus-Bolter 
correlation as: 

/!,.= 0 . 023 ^^ ( 4 . 9 ) 

dh 

The factor is dependent on the void and is intended to take into consideration the fact that 
with increasing wall superheating, the number of bubbles formed will increase. describes 
the heating surface fraction in direct contact with the sub-cooled liquid. 


^Kt) 


eff 

^he 


^ = l-7t S n./?2 

\e i = 1 


= 1 - nTzR. 


( 4 . 10 ) 


Where is the surface covered by bubbles. A;,, is the heated surface, n is the number of 
bubbles per unit area and nRh is the projection area of a bubble on the heated surface. 
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(1) Nnclcation cavity is activated 
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(2) Bubble is growing 



QcondcrE&tion 




A\ \ \ X>"\\ \ \ \ 's 
(4) Bubble after departmre 


Bubble is 
condensing 






ondcns&tion 



(5) Reconstruction of destroyed thermal layer 






(6) Bubble bcfoce collapse 

(a new bubble will grow) 


Figure 4.3: Schematic description of heat transfer mechanism by growth and 
condensation of steam bubbles in sub-cooled boiling regime. 


Although the bubbles have statistically different radii, the calculation has been earned out with 
average radius for the sake of simplicity. 
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Hainoun [1994] has proposed as; 


B 


KK 


n a 




for a < 


16a 


'aw 


(4.11) 


^KD for a > (4.12) 

where a^^,, is the void fraction at the point when the bubbles become detached from the heating 
surface(OSV).Experimental data shows that a„^^ may be about 5%-10%(Rogers et a/., 1987). 


Subcooling at OSV: The bubble will be detached from the wall only when the drag and 
buoyancy forces are greater than the holding force. The point at which the bubbles are first 
detached from the wall (this is OSV: Onset of Significant Void) is determined by the sub 
cooling of the liquid at this location. The sub cooling has been calculated using Saha and 
Zuber’s model [Saha and Zuber, 1974] as follows: 

= O.OOllq' 7^ Pe < 70000 (4. 13) 

\^t // 

/f 

=153.8 — ^ Pe> 70000 (4.14) 

Where Peclet number is defined as Pe = Re * Pr . 

Validity of the correlation; 

p = 0.1-13.8 MPa, m = 95-2760 kg.m'V, / = 0.28-1.9MWm-l 


Bubble Departure Radius: The radius at which the bubbles are detached from the wall after 
reaching their critical size and required subcooling is called the bubble departure radius. 
Bubble departure radius has been calculated using Unal’s semi empirical model [Unal, 1976]. 


R 


hd 


1.21x10“'' 



(4.15) 
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where 



d) = 



for 1/, >0.6W' 


= 1 for V,<0.61m5'"‘ 

_ r,-r, 

g^evap < ^ > 


(4.15a) 

(4.15b) 

(4.15c) 


The range of validity of the correlation: 

Pressure: 0.1-17.7 MPa, Heat flux: 0.47-10.64 MWm■^ Velocity: 0.08-9.15 ms'^ 
and Ar„„=3-86K 


4.3.3 Bubble Formation Rate 

The heat flux required for evaporation which is transferred from superheated boundary layer 
into the bubble can be calculated by: 

(iLp = « /n Vbd Ps Kvap (4-16) 

/„ is the detachment frequency of the nucleation center at which bubbles are formed and n is 
the number of bubble nucleation center per unit area of heating surface. The heat flux required 
to reconstruct the superheated thermal boundary layer is given by: 

- ^fn Qrecons (4-17) 
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Qrecons the quantity of heat per nucleation centre withdrawn from the superheated wall to 

reconstruct the thermal boundary layer. To eliminate the unknown product nf, a parameter E 
is being defined as: 


Er 


nfnQr 


Qgv _ , 

Qevap f 'I ^bd P g ^^evap 

If It 




tevap 


Hence, Qevap 


(4.18) 


Also from Meister’s consideration for thermal boundary layer [Meister, 1979], the parameter E 
has been calculated as: 




T.-T, 


^cvap - Evaporation parameter 


0.5 


(4.19) 


Finally, evaporation heat flux is given by: 




evap 




(4.20) 


Bubble formation rate: 

The bubble formation rate per unit area is given by: 


* ^evap 

* ~li 

evap 


(4.21a) 


Hence the bubble formation rate per unit volume can be calculated as: 

rr ^ 

r =r * *^ = r * 

* * A * /i„ 


(4.21b) 


evap 
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4.3.4 Bubble Condensation Rate 

The condensation of bubble in a sub cooled liquid is governed by two effects: (1) heat transfer 
at the phase interface, (2) inertia of the surrounding liquid. In the case of a large bubble formed 
with low sub cooling and low flow rate, condensation proceeds very slowly. Hence the inertia 
of liquid surrounding the bubble can be neglected. In this case condensation is mainly 
governed by heat transfer at the phase interface. On the other hand, small bubbles, resulting in 
the case of considerable sub-cooling and high flow rate, condense very rapidly. Owing to the 
inertia of condensation, the surrounding fluid cannot flow fast enough to fill the space vacated 
by the condensed bubble. Hence it cannot compensate for the resulting under-pressure. So, a 
large local pressure fluctuation may arise; this is termed as cavitations effect [Hamit, 1980]. 
Mayinger and Nordman [1976] have shown that the pressure fluctuations in the vicinity of 
condensed bubbles increases greatly above Jacob number 100. This is an indication that 
condensation is predominantly controlled by inertia above Jacob number 100. Chen [1986] has 
demonstrated by measurements on bubble condensation with the aid of holographic 
interferometry that a thermal boundary layer exists in the vicinity of condensed bubble, up to 
Jacob number of 60-80. Up to this boundary, the pressure fluctuations at the end of 
condensation are slight, which is an indication of heat transfer controlled condensation. Thus 
Jacob number can be used to differentiate different condensation regions. The Jacob number 
indicates the ratio between the energy that the liquid requires to reach the saturation state and 
the heat stored in the steam at the same volume. 


Ja = ) (4.22) 

0 1 

P ^evap 

(1) 7a < 80 : condensation is largely determined by heat transfer at the phase boundary. 

(2) 80<7a<100: transition region. Both the heat transfer and inertia effects are significant. 

(3) Ja > 100 : inertia effect is dominant. 
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Heat transfer controlled condensation 

Heat transfer controlled condensation rate has been calculated using Hainoun et aVs model 
[1996] as: 

3.6--^p° u, M/, Ja for Re^.^ < 10“ 

^bd 

r,. = 3.6 p° D; Nu, Ja for > 3x10“ 

‘^bd^h 

EeveRe^i, is the channel Reynolds’s number. 


(4.23) 


(4.24) 


For 10“ < Re^^ < 3x 10“ , F. is interpolated between the two regions. 

C, is the condensation parameter and equal to 0.16. Nu^ is given by Hewitt et al [1990] as : 


Nu, =0.185R4'>''“'^ 


(4.25) 


Rcg is the Bubble Reynolds’s number given by: Reg = 


YAm. 


(4.26) 


Nu^ = 0 . 228 Ref;iPr°-^A°-^^ 


(4.27) 


. has been introduced by Avdeev [1986] which can be found from the correlation: 

A,, =1 for a <5% , A^, = (l - for a > 5% (4.28) 


Inertia controlled condensation 

Inertia controlled condensation has been estimated by Hamit s correlation [1980] as follows. 
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r,. 




(4.29) 


is the condensation time [Hamit, 1980] which is given by: 


= 0.458c? 


bd 


kP‘j 


(4.30) 


Now the expressions for the evaporation rate in kgtn^s-' [equation (4.21b)] and 
condensation rates in kgm s [equations (4.23), (4.24) and (4.29)] have been calculated. 
Finally and F,. are put in the vapor mass balance equation [equation (4.2)] to find out the 
void fraction. 


4.4 Results and Discussions 


The comparison of simulated results with McMaster test data on axial void distribution has 
been given in Figure 4.4. The simulated results are matching quite satisfactorily with the test 
data. The length of heated section is 300 mm. Hence the void fraction keeps rising in between 
the length 100mm to 400mm as expected. As the rest of the section of the test pipe is unheated 
the void fraction is decreasing after the heated section. Figure 4.5 shows the channel pressure 
variation with respect to time at a low system pressure (1.542 bar). Figure 4.6 shows the 
channel pressure with respect to time at a high system pressure (72 bar). The average void 
distribution in the channel at low pressure and at high pressure is given in Figure 4.7 and 
Figure 4.8 respectively. 

The possible explanation for the void fluctuation at low pressure.may be described as 
follows. For the same quality, void fraction is larger at high pressure as compared to low 
pressure. This is because, the latent heat of vaporization decreases with increase of pressure. 
So, a slight change in quality causes a large change in void fraction at low pressure than at 
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high pressure. Hence, there is a greater chance that void will oscillate at lower pressure than at 
higher pressure. 

The reason for pressure fluctuation a, low system pressure may be attributed to the 
sudden collapse ol bubbles in sub cooled water above the heated section (riser) At high 

system pressure there is a less possibility of sub cooled boiling to occur. Hence the pressure 
remains constant at high system pressure. 



Figure 4.4: Comparison of simulated results with McMaster test data for axial void distribution. 
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Figure 4.6: Channel pressure variation with respect to time with initial pressure of 72bar. 
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Figure 4.7: Channel average void distribution with respect to time with initial pressure of 1.542 bar. 
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4.5 Closure 


A model for vapor generation and vapor condensation has been presented in this Chapter. The 
validity of the model has been checked successfully by validating the results with that of 
McMaster test data. The model is then used to simulate the geysering phenomenon at start up. 
In order to gain more confidence in the model, the code has been validated with the Apsara 
loop experimental test data, which is being done at Bhabha Atomic Research Centre (B ARC), 
Mumbai. This has been described in the next chapter. 
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Chapter 5 


Model VaHdation 

5.1 Apsara Experimental Loop 
5.1.1 Introduction 


Apsara reactor is a swimming pool type research reactor. It uses enriched uranium-aluminium 
alloy as its fuel material. It uses light water as moderator and coolant. Cadmium is used as 
control rod in this reactor. The core size of this reactor is 8.5 m (L) x 3.0 m (W) x 8.2 m (H). 
The maximum reactor power is 1 MW. 

Apsara loop experiment is being carried out at the Apsara reactor site. The neutron 
source required for the flow visualization is provided by the Apsara reactor. The objective of 
the experiment is to generate data for the validation of flow pattern transition criteria and 
developing correlations for flow pattern specific pressure drop. The experiment is being 
carried out jointly by Reactor Engineering Division (RED), High Pressure Physics Division 
(HPPD), and Solid State Physics Division (SSPD). HPPD has provided the expertise for the 
visualization of flow, recording of flow pattern and measurement of void fraction. SSPD has 
provided the support for neutron radiography. RED is responsible for the design, fabrication, 
safety and data acquisition. Installation at site and commissioning has been done by RED with 
on site help from SSPD and ROD. 
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5.1.2 Experimental Procedure 


First the total system is filled with water with the help of hand pump and vented out any 
stranded air inside the water. Some amount of water is drained out to maintain a certain level 
in the steam drum. Also pressurized nitrogen gas has to be supplied to bring the loop to the 
desired pressure level. The usual practice is to keep the pressure up to 4 bars higher than the 
desired system pressure for taking into account the leakage and dissolution of nitrogen with 
water. Once it is assured that the system is devoid of air bubble, the power is supplied to the 
system from heater and temperature is raised to the saturation temperature corresponding to 
the system pressure. Initially low power is supplied to the system. Later on as the temperature 
of the heated section as well as the fluid temperature stabilizes, power is increased. This raise 
in power is done in steps of 150W to 200W at a time. Maximum power to the heated test 
section depends on at what pressure the experiment is being carried out. After some time the 
system natural circulation starts, which becomes evident from the Ap readings. But two-phase 
will not be initiated as the actual pressure is still higher than the system pressure. To have two- 
phase natural circulation, some amount of pressure from the system is needed to be released. 
This is done through the vent valve provided at the top of the heat exchanger. Once the fluid 
temperature crosses the saturation temperature corresponding to the experimental pressure, it 
can be said that two-phase flow has been established. The experiment has been carried out for 
various tube sizes. The loop inventory corresponding to that tube size is given in the following 
Table 5.1. 


Test section(dia.) 

Inventory with condenser(liters) 

Inventory without condenser(liters) 

3/8” (7.035 mm) 

3.75 

2.072 

1/2” (10.21 mm) 

4.14 

2.460 

3/4” (15.74 mm) 

5.16 

3.490 

1” (19.86 mm) 

6.24 

4.560 


Table 5-1: Loop inventory corresponding to various pipe sizes 
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Maximum operating temperature and pressure of the system are 315 °C and 125 bar 
respectively. But in the present loop, the experiments were carried out at low pressure and 
temperature. Experimental set up is shown in the Figure 5.1. 
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TEST LOOP FOR FLOW PATTERN TRANSITION STUDIES 



SIMPLIFIED FLOWSHEET 


DESIGN PARAMETERS 
pressure : 125 bar 
Temperature : 315 °C 
Power : 10 kW 

Location : Apsara Reactor Hall 


Figure 5.1: The Apsara experimental loop. 


-Water Inlet 
-Water Outlet 
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5.2 Void Fraction Measurement Techniques 


Visualization of flow pattern and measurement of void fraction are very important to the study 
of two-phase flow. Observations of flow pattern are necessary to understand the physics of 
two-phase flow and to validate the flow pattern transition criterion. Conventionally optical 
method using glass tubes supplemented by X-ray, gamma rays and conductance probe 
methods are used for flow visualization. In the present work neutron radiography and 
conductance probe has been used. 


5.2.1 Neutron Radiography 


Most of the real life application in industry and particularly nuclear industry demands that 
flow pattern must be studied inside metallic pipe and that too under high pressure, temperature 
and heat flux. Conventional optical methods for flow visualization cannot be applied in such 
cases. X-ray or Gamma ray probe though easy to use, are not suitable for the simple reason 
that attenuation of these rays due to metallic wall of pipe is much more than that due to the 
water steam mixture. Also this technique is much less sensitive to the changes in the 
composition of the water-steam mixture. 

In the present neutron radiography thermal neutron has been used. It’s because thermal 
neutron can easily penetrate most of the metals used for pipes and have quite high sensitivity. 
Experimental set up of neutron radiography is shown in the Figure 5.2. The neutrons are 
absorbed by the water layers but passes through the voids. These unabsorbed neutrons are 
converted into light by the converter and it falls on the mirror where it gets reflected to another 
mirror. This is then viewed by the high sensitivity camera. Images are stored in image recorder 
and signals are passed to the data logger. 
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Camera controller 


Figure 5.2: Schematic of void fraction measurement by neutron radiography. 


5.2.2 Conductance Probe 


The schematic of conductance probe is shown in the Figure 5.3. It works on the principle of 
change of resistance with change in flow pattern. Working principle of the conductance probe 
is shown in the Figure 5.4 and Figure 5.5. Conductance probe is supplied with DC power 
supply, voltmeter and a resistance. When there is no void in the flow across the probe then the 
voltage measured is due to the resistance Ri. But when there is change in flow pattern i.e. 
voids in the flow, there is a large resistance across the probe than Ri. This means that the 
voltage across Ri is less than the former case. Conductance probe assembly is connected to the 
PC for data acquisition and processing. The void fraction is measured with the help of 
following formula. 
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Void fraction = 


Time elapsed between the changes of resistance 
Total time 


Ti 

Time 


Sealing assembly g g 



Figure S3: Schematic of Single point conductance probe. 


FLOW 



Figure 5.4: Electrical circuitry of a Single point conductance probe. 
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Figure 5.5: Voltage signal from the probe. 
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5.3 Results and Discussions 


The simulated results have been compared with the experimental data. The experimental as 
well as the simulated results have been reported below. 



Figure 5.6: Experimental single phase pressure drop in the horizontal leg of Apsara Loop. 


6000 
5000 
4000 
S 3000 

i 2000 
I 1000 
i „ 

-a 

g -1000 

I -2000 

-3000 



-4000 

-5000 


i I . . . . J 1 ^ 1 1 i 

250000 500000 750000 

Time in ms 


Figure 5.7: Simulated single phase pressure drop in the horizontal leg. 
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Figure 5.8: Variation in simulated single phase pressure drop with change in system pressure. 


Datc:02/06/03 for system pressure=sl9.22 bar 


•** vertical leg 
- horizontal leg 



Figure 5.9: Experimental Void fractions measured in Conductance probe. 
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Alpha-instanteneous 

Alpha-moving 

Alpha-weighted 

TEST-SECTION. after heater 
Power=3.637kw T-in=235, T-out=260, 
Pressure=51.8 bar 
Delta-p drop=30.4 mm of water. 
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Time frame no.[ 4 frame/sec] 


Figure 5.10: Experimental Void fraction measured by Neutron Radiography. 
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Figure 5.11: Simulated Void fraction in the two phase leg at system pressure 19.22 bar. 
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5.4 Code Tuning Parameters 

5.4.1 Numerical Stability 

The following parameters are essential in fine tuning the code. 

i) Upstream differencing parameter P : 

Convection Dominance: 


(a) A value of 0.2-0.3 for^ gives good results with LOCA conditions and at high 
pressure simulation. 

(b) = 1 gives full up streaming or forward differencing. 

P = 0 gives central differencing. 

The order of accuracy of central difference scheme is o(h^] whereas the order of accuracy of 
forward and backward difference scheme aieo{h). Therefore /3 = 0 will give a more accurate 
result; but this poses a problem in the stability of the numerical scheme. This is discussed in 
the following paragraph. 

= (5.1) 

^2<{> 

where a 2^ is thermal diffusivity and Pe number indicates the strength of convection and 

diffusion in two phase flow. Hence, the condition Pe » 1 indicates that convection is more 
dominant. The qualitative temperature distributions T vsx (with Pe as parameter) are 
indicated in Figure 5.14. 
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Figure 5-14: The fluid temperature profile as a function of Peclet number. 


From the above figure it can be inferred that when convection is more dominant then upstream 
differencing scheme is more realistic as it uses backward difference for the convection term. 

In the present study, the convection term is more dominant as Pe»\. Hence ^ =0.8 -0.9 is 
suggested; since ^ more close to 1 means proximity to upstream differencing. 



Large Grid 


Small Grid 


Figure 5.15: Grid Size 


When using smaller grid size upwinding is more realistic as RHS value will clearly depend 
upon LHS. On the other hand, central differencing is suggested when using larger grid size. As 
the difference in values at LHS and RHS is large, an average of both the values will fit as the 
value of the entire cell. 
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ii) Time Step Ar 

Apart from calculating the effective Courant number At can be utilized to choose a proper 
method out of Pure-implicit, Implicit-explicit (semi implicit) and explicit schemes. Courant 
number is defined as At/Ay . 


T? 



Explicit 


Figure 5.16: Physical representation of various numerical schemes. 

Hence if At is large, then pure implicit is more realistic whereas for small At semi-implicit is 
more realistic. 

iii) Stability Control Parameter. 

It is used to set the effective Courant number for stability of the numerical scheme as 

^ 1 _1 (5.2) 

Ay STABC 4 

The purpose of viscosity in viscous stress term is to damp the high frequency oscillations that 
can not be resolved or would lead to numerical instability. But the magnitude of Oj (Specified 
Kinematic viscosity coefficient) must not be so large as to appreciably damp the larger scale 
flow instabilities of physical interest. 
iV ) At 

Let, ^ = f 

Ay 
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V)^ is calculated using Hirt’s[C.W. Hirt, 1968] results as: 


u, « 


A}- 

"J 'max 


2 


(^mLofAj Ay At 


Ay At 


Ay 


I (^J^At 
2 Ay 


Ay 


1 (^m) A/ 

|_ ^ m /max 

2 Ay 


M. 

2At 


(Ayf 


2At 


To avoid diffusion instability, we also require 


(5.3) 


4t)jAr < (Ay)* with 

P2(!> 

This is automatically satisfied provided / < 1/2 . 

/ = l/4 provides a better solution, which simultaneously ensures sufficiently small fluid 
motion per cycle and adequate dissipation, while precluding diffusion instability and excessive 
damping. 


iv) Pressure Convergence Rate Parameter. 

This is the relaxation parameter (OMG) for the calculation of Ap in the pressure iteration. 
This variable affects the rate of convergence. A value of OMG = 1.7 often improves the rate 
of convergence for low Mach number flows. For high Mach number flows (Supersonic), 
OMG = 1.0 or near to 1 is suggested. Two phase flow under natural convection is at Mach 
number less than unity. 
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v) Convergence Criteria. 

While working with low pressure more stringent (eslO"^) convergence criteria can be set. 
Also while working with double precision, a strict convergence criteria (e = 10”*) is expected. 
For single precision, generally the convergence criterion is of the order of(e = 10”*). However, 
more strict criteria lead to more CPU time for mnning the code. 

A sensitivity analysis has been carried out in double precision to study the effect of £ 
on the results and on the number of iterations to converge. 



Figure 5.17; Pressure Variation with £ as a parameter. 
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Figure 5.18: Variation of Number of Iterations with 8 . 


From the sensitivity analysis it is concluded that applying a stricter convergence criterion does 
not refine the results. On the contrary it increases the number of iterations for convergence; 
thereby increasing the CPU time. 

vi) Subcycle Parameter. 

This specifies the desired number of sub cycles in the code. A large network system often 
contain low speed flow (Low Mach Number) with slowly varying properties in one region and 
high speed flow or flow that requires a finely detailed description in another region. Hence a 
provision is made in the present work that allows the use of different time steps for integration 
in each component and junction cell in the network! It’s because the time step can be 
significantly different in the various components; that is one component may have to be 
integrated several time steps to keep pace with one time integration in another component. 

This subcycling feature enables accurate time integration. 


95 




Model Validation 


5.4.2 Heat Transfer 

i) Heat Flux and Wall Temperature. 

Q/tux - calculated directly from the power supplied to the test section. As the power 

that is supplied has to be distributed uniformly in the test section volume, the heat flux per unit 
volume is given by 

_ Power 

^^flux 2 

%*r. * Height of the Test sec tion 


(b) The second option may be to calculate thej2^„ from the given wall temperature of 
the test section. 

= ft, AAr, (5.4) 

Dittus-Bolter correlation gives 

h. = 0.023-^ PrJ-'' where Re = and Pr = . (5.5) 

4, k„ 

(5.5a) 

J_ = ^ + iz£ (5.5b) 

A,. = Ttd,. * Height of the Test sec tion 
AT.=T. -T 

I wi m 


ii) ETEM 

ETEM=1; Liquid and Vapor phases are maintained at equal temperature. 

ETEM=0; The vapor phase is maintained equal to saturation temperature at the local pressure. 
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Figure 5.19: Phase equilibrium diagram on T-S coordinate. 


r, = Tg essentially means it can be at state 1 or 2. 

T„ =T means it can be at 3 or 4. 

g A 

When inventory is very large and void fraction is less, then the large heat content of the 
liquid phase keeps the liquid temperature nearly invariant. For such case T, = or 

makes little difference. For low flow and low pressure condition with sufficiently large void 
(a > 0 . 5 ) the two phases move with a very small speed. Hence the interaction time between 
the two phases is large. Therefore the possibility of attaining equilibrium is bright. On the 
other hand, in the case of high pressure, high void with turbulent flow condition (High 
Reynolds number), is a highly non-equilibrium state. Here the assumption of saturated vapor at 
local pressure will fit well. 
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iii) Number of Bubbles per unit Volume 



Figure 5.20: Vapor Bubble Nucleus 


If /?^=Radius of curvature of the nucleus surface, then surface tension creates a pressure 
within the nucleus such as 


P, =Pi+- 


Clausius-Clapeyron equation gives 


dT J I 1 


pI Pi 


Assuming that the steam behaves like an ideal gas, -4- « -V- Hence is neglected. 

P? P° Pi 

Therefore = (5.8) 

dT T 

2.0 T 

Integrating and using equation ( 1 ) we have = — r — (5 .9) 

Pgh„,p -T^ 

is also known as active nucleus radius; that is only after reaching R^, nucleus will be 
activated and can lead to a bubble. 

From equation (2) it is clear that with increasing wall superheat ( t^ th® active 

bubble radius {R^) will decrease. That means all nuclei of higher radii than R^ will be 
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activated. Thus the chances of formation of large nucleating site have enhanced. This is also 
confirmed by the following fact. 

=l-n,nR- 

where is the heating surface fraction which is bubble free and in direct contact with sub- 
cooled liquid. TtR^ is the projection area of bubble on the heating surface. will decrease as 

the radius of active nucleus R^ decreases. Therefore, for a constant A,,, with decreasing 
will increase. 


iv) p® AS A FUNCTION OF TEMPERATURE 



0 10 20 30 40 50 60 70 80 90 100 


Figure 5.21: Variation of density with temperature and the residuals of quadratic fit. 
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The variation of liquid density with temperature has been plotted. Then a quadratic curve 
fitting has been performed to obtain the equation of dependence of density with temperature. 
The equation so obtained is found to be, 

p? =-0.0036*r^ -0.076*7+1000.0 (5.10) 


5.4.3 Hydrodynamics 

i) Relative Velocity Parameter 

DFVEL is a program control parameter that determines whether the relative velocity is to be 
calculated or not. A value of DFVEL =1 calculates the drift velocity whereas DFVEL = 0 
bypass the drift velocity subroutine. 

Hence, a value of DFVEL = 0 and ETEM = 1 will convert the Drift Flux Model to a 
HEM (Homogeneous Equilibrium Model). 

ii) Turbulence Parameter p, 



= Mass averaged mixture velocity 
= Average relative velocity 
= Parameter accounting for turbulent fluctuation 
Generally p, has a value of 0. 1 or less because large turbulent velocity fluctuations often have 
magnitudes as large as 10% of the mean velocities. Nevertheless, the best value of P, must be 

determined by comparison with experimental data. 

A sensitivity analysis of P, has been performed with respect to the mixture velocity. 

From the analysis it can be concluded that P, has a very minute effect on the mixture velocity 
at high system pressure. 
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Figure 5.22: Variation of Mixtu 



Figure 5.23: Enlargec . . 




Conclusions and Scope for Future Work 


Chapter 6 

Conclusions and Scope for Future Work 

6.1 Conclusion 


A four equations Drift Flux model has been used for numerical simulation in the present work. 
The focus of the research is to bench-mark the geysering phenomenon. In the first stage of the 
present work the code has been tested to simulate small break LOCA condition. In the second 
stage, a vapor generation and bubble condensation model has been incorporated in the code. 
The model has been then validated with McMaster experimental test results on axial void 
fraction distribution. This model has been used to simulate geysering phenomenon in the same 
McMaster experimental set up. Finally the usefulness of the code has been studied while trying 
to validate the Apsara Loop experimental test data for void fraction at the two phase leg, 
pressure drop in the single phase leg and temperature of the liquid in the loop. 

6.2 Scope for Future Work 


The present work has opened up new areas of exploration. Following challenges are 
recommended for future researchers. 
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1 . The simulation of Apsara Loop has been done by considering two sub-systems. One 
consisting of single phase leg and the other having two phase flow. The simulation of 
Apsara Loop as a whole is desirable to obtain more realistic solution. 

2. A non-uniform initialization of all primitive variables (Pressure, Temperature, Velocity 
and Void fraction is expected. The code should have the provision for taking 

(a) initial condition as given and 

(b) initial condition as output of previous simulation. 

3. Temperature dependent initialization of microscopic liquid density (p° ) should be 
incorporated. 

4. The coalescence of bubbles has not been included. Therefore a model for the 
coalescence of bubble should be incorporated. 

5. The model should be extended to two dimensional and three dimensional domains for 
better visualization. 

6. Better estimation of momentum exchange function has to be made taking into account 
Basset, Virtual mass and Faxen forces. 
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Appendix 

A.l Void Fraction 

The macroscopic densities for vapor and liquid are related to the microscopic quantities as, 

p^, =apj, P;=(l-a)p? 

p„, =apj, +(l-a)p° 

= p^, +p?-ap? 
ap°=p^+p°-p„, 

where p° = Microscopic vapor density or density of the pure vapor 
p° = Microscopic liquid density or density of the pure liquid 


A.2 Vapor Continuity Equation 


The vapor mass balance equation is given by 


9pp 1 3 




dt A dy 


PpP/ 






Continuity equation is discretized at the cell centre j 
Hence discretizing Equation (A.l) at j 


(A.1) 
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Also since p,„K„ +p,K, and7, =V;, -V,, 

Now p^V„,+- 
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^jh-v 
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p*‘'«(p*+p/i 


Pm '' Pm 


But p„, =apj, +(l-a)p{’ andp^, =ap;, p, =(l-a)pf 
Thatis, p„, =p^, +p, 


P.P/ 


P/f'Ki.’P'i 


Hence p^V,„ +^V, = ' '- = p^V, 


Now — -^/l] 
A dy 


f 


Pt-'^'m + 


P.P/ 


Pm 




'g - m ■ *'r 


=tI:'^(p.''J 


J 


A dy 


Discretizing Equation (A.2d) at J 


11 . 

Ady 




= -i-A 

. = A,a3. 


U(pA>)l 


Applying central differencing scheme the RHS becomes 


I J 2 


Aj Ayj 

This is the form of equation reported in Equation (2.25). 


(A.2a) 


(A.2b) 


(A.2c) 


(A.2d) 


(A.3) 
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A.3 Mixture Energy Equation 


The Mixture energy balance equation is given as, 


I ^ ^ /(I 

3; a'^^' 


0 IV +' ^ 


P.P/ 




V. 

=-^i-A 

y,n + 

r I 1 ) 

/ r 

A dy 

Pm 

[p/ Pt°J 


+ Kv;+W,^+Q 


The Mixture mass conservation equation is 

^P,n , 1 


■ + ■ 


■ = 0 


dt A dy 
Multiplying Equation (A.5) by /„, gives 

, ap,. , L 3(Ap„v.,) „ 

" 3r A dy 

Subtracting Equation (A.6) from Equation {A.4) yields 


3P., , „ /„ 3(Ap,V,J Ap„V„ SI, _1 ^ , 

"ir’^P’lT'^T 3), “T” 3y A3/' 


j ''t-'m 


_{/ ^Pm , An 


dt A dy 
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P.P/ 


1 1 
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(in-l.y. 


p,° p/° 
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J J 


+ KVr +Kl 


dl. 


P„,^ + P,n^,nM- + T-^l 


M + IA, 

dy A dy 


PjfP; 




= - £ Aa 

A dy 


T. PsP' 


P« 


1 I 

0 ^ ( 


Pt 


+ KVj'+W,i,+Q 


Dividing throughout by p,„ in Equation (A.8) gives, 


(A.4) 


(A.5) 


(A.6) 


+ Q (A.7) 

■N ■ 

/ 

(A.8) 


111 



Appendix 


9/ 3v p,„Aav i 





, P.P/ 

P,« 




Pm Pm Pm 


Mixture energy equation has to be discretized at the cell centre j . 
Hence discretizing Equation (A.9) at J will give, 




(A.9) 

(A. 10a) 

(A, 10b) 



(A. 10c) 


(A.10d) 


(A.10e) 
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—kv/+—w,,,+—q 

P m Pin P in 



(A.10f) 
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The above Equations (A.9) to (A. lOf) will give the updated value of mixture energy (/^ 
which can be used as a consistency check for the Equation (2.35). 


A.4 Two Dimensional Equations for Junction Cell 

U = Velocity in the X - direction 
V =Velocity in the Y -direction 


Total Mass Balance Equation 


^Pm I 1 

dt A 


• a(Ap,„t/„,) ^ 3(Ap,.vJ ' 

dx dy 


(A.11) 


Vapor Mass Balance Equation 


dt A 


dx‘ 




+— Ai 

oy 


p y + 

r> '■ 


= r -r 

s <■ 


(A.12) 


Mixture Momentum Equation in X-direction 

^ n n. f 


3p,U„. , 1 

at A 


dx 


2 . P«P/ rr 2 
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Pm 
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dp f 
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(A.13) 


Mixture Momentum Equation in Y-direction 
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dp 

f v is 


(A.14) 


Mixture Energy Balance Equation 


^9,Jm 

dt 


+ ■ 


dx 
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P,P/f 1 







v;„ + 


p«p/f 1 




Vr 


+ k[u/+V;\+W,^+Q (A.i5) 


A.5 Variables for Constructing the Junction and BC 

KBOT (I): Number of junction connected to the bottom of pipe I. A zero value indicates an 
isolated end. 

KTOP (I): Number of junction connected to the top of pipe I. A zero value indicates an 
isolated end. 

LANG (I): Parameter that indicates the orientation of component I. Refer Figure 3.4 and 
Figure 3.5. 

LBOT (I): Parameter that indicates the type of boundary conditions to be applied at the bottom 
of component I when it is an isolated end. 

LTOP (I): Parameter that indicates the type of boundary conditions to be applied at the top 
of component I when it is an isolated end. 

MBOT (I): Parameter that specifies the boundary data set for the bottom end of component I. 
MTOP (I): Parameter that specifies the boundary data set for the top end of component I. 
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1. INTRODUCTION 

The current Light Water Reactors (LWR) achieves shut 
down through active safety systems. Passive safety 
systems have been proposed for advanced designs to 
enhance the reliability of safety functions. The 
advanced LWR would incorporate a number of passive 
safety features in its design. One of them is to adopt 
natural circulation core cooling during start-up, power 
raising, rated power conditions and accidental 
conditions. This concept is to eliminate the recirculation 
pumps which are normally present in conventional 
forced circulations BWRs. This is due to the reason that 
the forced circulation - loop has the disadvantage of 
using a pump, which is costly. Again if by any reason 
the pump breaks down then the dissipation of fission 
heat is hampered. This results in tremendous 
accumulation of heat and may cause core melting in the 
reactor. Thereby it demands backup safety measures 
which will add to the cost of reactor. 

However natural circulation systems require power to 
initiate the circulation through void generation. This 
i;nean the natural circulation reactor would be heated by 
fission energy from the startup under low temperature 
and low pressure condition. Thermal-hydraulic 


instabilities have been reported under low pressure 
conditions (Chiang et al, 1994). If thermal hydraulic 
instabilities were to occur at startup then the reactor 
would not potentially continue operation during 
power up because the void fraction fluctuation in the 
reactor core would oscillate the reactivity. Therefore 
it is necessary to investigate and understand properly 
the thermal hydraulic instabilities during start up. 

Aritomi et al. (1992, 1993) and Chiang et al. (1992, 
1994) have conducted extensive research in Llic aacii 
of geysering under natural circulation. Masuhara et 
al. (1993) have also performed small scale 
experiment to demonstrate this phenomenon. These 
experiments illustrated that the geysering mode 
oscillation would occur at low pressures and low 
flow conditions, Aritomi et al, (1992) had explained 
the driving mechanism of geysering as follows: 
When voids are generated in a heated channel, a large 
slug of bubbles forms, which grows due to decrease 
in hydrostatic pressure head as it moves towards the 
exit. The vapor then mixes with the liquid is the 
subcoolcd riser or upper Plenum and is condensed 
there. Due to bubble collapse and subsequent 
decrease in pressure, the subcooled liquid 
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inters the channels and restores the non-boiling 
iditions. This process repeats periodically causing 
w oscillations. Hence it is evident that the bubble 
maiions, growth and collapse phenomenon are of 
;>onance to geysering instability. 


3. GOVERNING EQUATIONS 

The basic one-dimensional, four equation drift flux 
model for two phase flows consists of the following 
conservation equations 


rlier attempts to model startup instabilities (Aritomi 
iL, i992)and (Paniagua et aL, 1996) indicate that, to 
diet the possible startup instabilities correctly, it is 
portant to accurately predicts the vapor generation 

In the present code the effect of the bubble 
malion a.s well as the condensation rate has been 
isidered. Moreover, most of the models are 
dsaged for system with high pressures (greater tlian 
bars) and thus are unsuitable for the simulation of 
^sering owing to the great influence of pressure on 
void content in the subcooled boiling regime, 
ally in the existing models there is no conclusive 
I physically well-defined description of the mass 
isfer rate between the vapor and liquid phase, an 
ect which is of great significance for the 
ierstanding of phenomena in the subcooled boiling 
ime. In the present code a four equation drift flux 
del has been used which is numerically more stable 
n the five equation and six equation model. 

2: FORMULATION AND SIMULATION 
SCHEME 

i following code is based on the following thermal 
Iraulic modeling features: 

bur fundamental balance equations - one liquid 
5s balance, one vapor mass balance equation, one 
lure momentum equation and one mixture energy 
mce equation. 

econd order finite difference formulation. 

>rift flux model for phasic velocities 

appropriate vapor generation and vapor condensation 

lodel. 

h component has a one-dimensional representation 
1 a variable cross-secjtional area. The equations are 
'ed by a partially ^plicit method that can use 
erent time steps in different components. The 
iponents are discretized using staggered mess 
ngements. The momentum equation is advanced 
licitly that is, explicit updating of velocity has been 
e. The pressure has been updated implicitly. When 
transients are initiated, the system pressure changes 

assumed to be instantaneous and uniform 
ughout the individual computational cells. This 
is to a quasi-steady pressure distribution throughout 
system. The objective of our code is to determine 
void fraction in the subcooled regime with proper 
^deration of bubble formation rate and bubble 
lensation rate. 


Liquid mass balance: 


dt 




u 


he 


c A g 


( 1 ) 


Vapor mass balance: 


^{ap VALp V 1=- 
aA &) dyV g g) 


u. 

■ Tc +— 

A g 


( 2 ) 


Mixture momentum balance: 


— (p V ]+—[ccp V ^ +{\-a)py^'\ 
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+ Pgy-^fvis 


Mixture energy balance: 



= (4) 

Here is the distributed losses (pipe wall friction 
and local losses due to sudden change in area), U 

is the heated perimeter, 1 ^ is the mixture specific 
internal energy: 

Pm^rn is the energy 

dissipation term, is the relative velocity between 
phases: and Q is the heat source term 

in kW. 


4. MODELLING OF VOID FORMATION IN 
THE SUBCOOLED BOILING REGIME 

In Ihc subcoolcd boiling regime, boiling occurs at the 
liquid cooled heating surfaces due to high heat flow 
densities, although the fluid has on average not yet 
reached the saturation temperature associated with the 
system pressure. In accordance with a suggestion by 
Griffith (Griffith et al, 1958) four zones of heat 
transport and flow activity can be differentiated along a 
channel axis until saturation boiling has been reached. 
Zone I is referred to as single phase heat transfer zone, 
Zone II is onset of nucleate boiling (ONB) zone. Zone 
III is onset of significant Void (OSV) and Zone IV is 
the saturated core flow region. The equation of the 
following effect must be formulated in order to 
determine the void content in the subcooled boiling 
regime: (1) on set of nucleate boiling (ONB). (2) bubble 
formation and bubble growth (bubble generation rate) 
at the heating surface (3)onset of significant void 
(OSV) and bubble departure diameter (4)bubble 
condensation in the subcooled core flow. 

4.1 Onset of Nucleate boiling (ONB) 

It describes the point at which first boiling nuclei are 
activated at the boiling surface. A correlation by 
Bcrgels and Rohsenow (Bergels et al, 1981) has been 
used to determine the wall superheating at which ONB 
will be activated. 
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The heat flux associated with the single phase 
convection from the heating surface is given as: 

9l<t) = ) 0) 


Where Tw and Ti are the wall and liquid temperature 
respectively and /ij^ is the heat transfer coefficient of 
the single phase liquid-flow, which can be determined 
from the Dittus-B enter correlation as: 




0.023-^ Re 
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The factor i?j<|,is dependent on the void and is 
intended to take into consideration the fact that with 
increasing with the wall superheating, the number of 

bubbles formed will increase. describes the 
heating surface fraction in direct contact with the 
subcooled liquid. 
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The range of validity of the correlation is: channel 
diameter 2.4 and 4.6mm; flow velocity: 3-7ms‘ 
‘.Pressure: 1-36 bar. 

4.2 Onset of Significant Void (OSV) 


Where is the surface covered by bubbles, is 
the heated surface n is the number of bubbles per unit 

area ^nd ttR^^ is the projection area of a bubble on 
the heated surface, Although the bubbles have 
statistically different radii, the calculation has been 

carried out with average radius for the sake of 
simplicity . 

Hainoun(l994) has proposed as: 


Assuming that the subcooled boiling occurs, the total 
heat Q supplied to the fluid from the heating surface 
per unit area can be arbitrarily partitioned into the 
following fractions: i 

2 = 2,o+Gav forex^-f- 

Heat required to generate void is given by. 
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whei:e is the void fraction at the point when the 
bubbles become detached from the heating 
surface(OSV). Experimental data shows that may 
be about 5%-l()%(Rogers et r//..1987) 

Siibcooling at OSV: The bubble will be detached from 
the wall only when the drag and buoyancy forces are 
greater than tiie holding force* The point at which the 
bubbles are first detached from the wall (this is OSV : 
Onset of Significant Void) is detennined by the 
subcooling of the liquid at this location. The subcooling 
has been calculated using Saha and Zuber’s model 
(Saha and Zuber, 1974) as follows: 


4,3 Bubble fonnatlon rate 

The heat llux required for evaporation which is 
transferred from superheated boundary layer into the 
bubble can be calculated by: 

q (15) 

/is the detachment frequency of the nucleation 
center at which bubbles are formed, n is the number 
of bubble nucleation center per unit area of heating 
surface. Tho heat flux required to reconstruct the 
superheated thermal boundary layer is given by: 


= 0.0022^'-^ Pe < 70000 (12) 
k, 

Jf 

Ar, =153.8 — ^ Pe> 70000 (13) 

ViP.Cp, 

where Peclet number can be found by the correlation 
Pe=Re 

Validity of the correlation: 

P=0.1-13.8 MPa. ni=95-2760 kg.m'V*. 

=0.28-1. 9MWm-^ 

Bubble ’Departure Radius: The radius at which the 
bubbles are detached from the wall after reaching their 
critical size and required subcooling is called the bubble 
departure radius. Bubble departure radius has been 
calculated using Unafs semi empirical 
model(Unal,1976) 
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The range of validity of the correlation: 
Pressure: 0.1-17.7 MPa 
Heat flux: 0.47-10.64 MWm'^ 

Velocity: 0.08-9.15 ins'* and 

Liquid sub-cooling: AT^^^ =3-86 K 


(Zo-V = ^fO-rccon. 


(16) 


Qrccons quantity of heat per nucleation centre 

withdrawn from the superheated wall to reconstruct 
the thennal boundary layer. To eliminate the 
unknown product nf, a parameter E is being defined 
as: 
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Also from Meistcr consideration for thermal 
boundary layer (Meisier, 1979), the parameter E can 
be calculated as: 


T = 2C 
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Finally, evaporation heat flux is given by: 
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Bubble formation rate: 
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(20) 



4.4 Bubble condensation rate 


a 


The condensation of bubble in a subcooled liquid is 
governed by two effects: (1) heat transfer at the phase 
interface; (2) inertia of the surrounding liquid. In tlic 
case of a large bubble formed with low subcooling and 
low flow rate, condensation proceeds very slowly. 
Hence the inertia of liquid surrounding the bubble can 
be neglected. In this case condensation is mainly 
governed by heat transfer at the phase interface. On the 
other hand, small bubbles, resulting in the case of 
considerable sub-cooling and high flow rate, condense 
very rapidly. Owing to the inertia of condensation, the 
surrounding fluid cannot flow fast enough to fill the 
space vacated by the condensed bubble. Hence it cannot 
compensate for the resulting underpressure. So, a large 
local pressure fluctuation may arise; this is termed as 
cavitation effect (Hamit, 1980). Mayinger and Nordman 
(1976) have showed that the pressure fluctuations in the 
vicinity of condensed bubbles increases greatly above 
Jacob number 100. This is an indication that 
condensation is predominantly controlled by inertia 
above Jacob number 100. Chen (1986) has 
demonstrated by measurements on bubble condensation 
with the aid of holographic interferometry that a 
thermal boundary layer exists in the vicinity of 
condensed bubble up to Jacob number of 60-80. Up to 
this boundary the pressure fluctuations at the end of 
condensation are slight which is an indication of heat 
transfer controlled condensation. Thus Jacob number 
can be used to differentiate different condensation 
regions. The Jacob number indicates the ratio between 
the energy that the liquid requires to reach the 
saturation state and the heat stored in the steam at the 
same volume. 
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for >3x10'^ 
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For 10"* < < 3x10"^, is interpolated 

between the two regions. 

where is the condensation parameter and equal to 
0.16. Nu^ has been given by Hewitt et al (1990) as : 

Nu, = 0.185 (24) 


Re^ is the Bubble Reynold’s number given by: 

r, 

ReB=-!-^ (25) 

Nu^ = 0.228 Ref Pr®-^ Av°“ (26) 

Here Re^ is the channel Reynold’s number. Av has 

been introduced by Avdeev (1986) which can be 
found from the correlation: 

Av=l for a <5% 

Av = (l — for cc > 5% 

(27) 

Inertia controlled condensation 

Inertia controlled condensation has been estimated by 

Hamit’s correlation (1980) as follows: 


r Pl^pjb'soi ^l) 

Ja = (21) 

p h 

r g evap 

(1) Ja < 80 : condensation is largely determined by 
heat transfer at the phase boundary. 

(2) 80 < Ja <100: transition region. Both the heat 
transfer and inertia effects are significant. 

(3) Ja > 100 : inertia effect is dominant. 

Heat transfer controlled condensation 

Heat transfer controlled condensation rate has been 

calculated using Hainoun et aVs model (1996) as: 

= C^3,6—^ pgViNu^Ja 
^bd 

for Re^^ < 10^ 
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(28) 


T^is the condensation time (Hamit, 1980) which is 
given by: 


r = 0.458J 
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Now the expressions for the evaporation rate in 


kgiif^S~^ [equation (20)] and condensation rates 

in kgrn '^ [equations (22), (23) and (28)| 

have been calculated. Finally F^ and F^. are put in 

the vapor mass balance equation [equation (2)] to 
find out the void fraction. The simulation has been 


( 22 ) 
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verified by comparing it with tlie McMaster 
dmental data on axial void distribution. 

5. RESULTS 



[ Comparison of simulated results with McMaster 
test data for axial void distribution in sub-cooled 
boiling regime 



I Channel Pressure Variation with respect to time 
vith initial pressure of 1 .542 bar 
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Fig. 3 Channel Pressure Variation with respect to time 
with initial pressure of 72 bar 



Fig. 4 Channel Average Void Distribution with respect 
to time witli initial pressure of 1.542 bar 







Fig. 5 Channel Average Void Distribution with respect 
to time with initial pressure of 72 bar 

6. CONCLUSION 

A computer code has been developed for the calculation 
of axial void fraction in the subcooled boiling regime 
using a four equation drift flux model with proper 
consideration of bubble formation and bubble 
condensation rates. The results obtained are in 
agreement with the McMaster experimental data on 
axial void distribution, as shown in figure 1. In figure 2, 
the variation of channel pressure with respect to time, 
with initial pressure 1.542 bar, is shown. The figure 
shows that at a low initial pressure during startup, 
boiling instabilities are prominent, which is a clear 
indication of geysering. The geysering period can also 
be determined from this figure, which is found to be 
about 20 ms for the case studied. Figure 3 shows that at 
a high initial pressure during startup, boiling 
instabilities are absent. The void fraction variation 
shown in figures 4 and 5 also confirm the above 
conclusion. 
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